August 2024: Plots based on Model Selection for Environmental Parameters in SGD (this script is used for majority of results)
Hannah Merges
2025-05-14
Loading libraries
Loading Data
Reviewed May 15th, 2025:
Plots for Gross Photosynthesis - Cabral:
- Cabral is significant for silicate at linear (0.04) and nonlinear (0.02)
- Cabral is only significant for SGD Dils in the nonlinear model (0.02)
- With removing an outlier point in the last dilution, Cabral is significant for GP both linear (0.015) and nonlinear (0.015)
# mutate(ammonia_umolL = ifelse(ammonia_umolL = "<0.02", 0.01, ammonia_umolL)) %>%
###################
## Silicate
####################
### confirm linear model
GP_Cabral_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
#filter(site=="Cabral") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Cabral_silicate_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(silicate_umolL + 1) 0.060043 0.060043 1 110.84 3.2591 0.07374 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### confirm nonlinear model
GP_Cabral_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(silicate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
# filter(site=="Cabral") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Cabral_silicate_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(silicate_umolL + 1), 2) 0.20901 0.10451 2 110.86 6.1685 0.002883
##
## poly(log(silicate_umolL + 1), 2) **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### both are significant so can add either trend line
GP_Cabral_Silicate <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=silicate_umolL,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~log(x)", color="black") +
# geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
# geom_hline(yintercept = 0) +
labs(x = "log silicate",
y = "Rate (umol O2 cm2 hr-1)",
title = "GP Rates for Cabral with Silicate") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
GP_Cabral_Silicate
###################
## SGD dils
####################
### confirm a linear model
GP_Cabral_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Cabral_SGD_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(SGD_number + 1) 0.045966 0.045966 1 46.338 3.266 0.07722 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### confirm a nonlinear model
GP_Cabral_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Cabral_SGD_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.10382 0.051908 2 45.272 3.9695 0.0258 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1GP_Cabral_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=SGD_number+1,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
# geom_hline(yintercept = 0) +
labs(x = "log SGD Dils",
y = "Rate (umol O2 cm2 hr-1)",
title = "GP Rates for Cabral with SGD Dils") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
GP_Cabral_SGDDils
###################
## SGD dils WITHOUT POTENTIAL OUTLIERS
####################
GP_Cabral_SGDDils_nooutlier <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(sample_ID!= "Cabral_Col6_Dil9_Light") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=SGD_number+1,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9)) +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
#geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
geom_smooth(method="lm", formula = "y~log(x)", color="black") +
# geom_hline(yintercept = 0) +
labs(x = "log SGD Dils",
y = "Rate (umol O2 cm2 hr-1)",
title = "GP Rates for Cabral with SGD Dils") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
GP_Cabral_SGDDils_nooutlier
### confirm a linear model
GP_Cabral_SGD_linearmodel_nooutlier <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(sample_ID!= "Cabral_Col6_Dil9_Light") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Cabral_SGD_linearmodel_nooutlier)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(SGD_number + 1) 0.081847 0.081847 1 50 6.2821 0.01549 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### confirm a nonlinear model
GP_Cabral_SGD_nonlinearmodel_nooutlier <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="GP") %>%
filter(sample_ID!= "Cabral_Col6_Dil9_Light") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Cabral_SGD_nonlinearmodel_nooutlier)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.11464 0.057321 2 49 4.5402 0.01553 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Plots for Gross Photosynthesis - Varari
- Varari is ONLY significant for silicate in nonlinear model (0.02)
- Varari has no significant relationship with SGD and GP
###################
## Silicate
####################
### confirm linear model
GP_Varari_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Varari_silicate_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(silicate_umolL + 1) 0.015463 0.015463 1 62.55 0.7003 0.4059### confirm nonlinear model
GP_Varari_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(silicate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Varari_silicate_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(silicate_umolL + 1), 2) 0.15527 0.077633 2 61.577 3.925 0.02487
##
## poly(log(silicate_umolL + 1), 2) *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#### ONLY significant in nonlinear model
GP_Varari_Silicate <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=silicate_umolL,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
# geom_smooth(method="lm", formula = "y~log(x)", color="black") +
geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
# geom_hline(yintercept = 0) +
labs(x = "log silicate",
y = "Rate (umol O2 cm2 hr-1)",
title = "GP Rates for Varari with Silicate") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
GP_Varari_Silicate
###################
## SGD dils
####################
### confirm a linear model
GP_Varari_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Varari_SGD_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(SGD_number + 1) 0.010297 0.010297 1 60.065 0.4631 0.4988### confirm a nonlinear model
GP_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Varari_SGD_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.089681 0.044841 2 59.102 2.1119 0.13GP_Varari_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=SGD_number+1,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
# geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
# geom_hline(yintercept = 0) +
labs(x = "log SGD Dils",
y = "Rate (umol O2 cm2 hr-1)",
title = "GP Rates for Varari with SGD Dils") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
GP_Varari_SGDDils
Plots for Respiration - Cabral:
- There is no significance at Cabral for any environmental parameter or broad SGD dilutions
###################
## SGD dils
####################
### confirm a linear model
R_Cabral_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="R") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(R_Cabral_SGD_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(SGD_number + 1) 9.037e-08 9.037e-08 1 46.085 0 0.9962### confirm a nonlinear model
R_Cabral_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="R") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(R_Cabral_SGD_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.011554 0.0057771 2 45.076 1.5607 0.2211R_Cabral_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="R") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=SGD_number+1,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
# geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
# geom_hline(yintercept = 0) +
labs(x = "log SGD Dils",
y = "Rate (umol O2 cm2 hr-1)",
title = "R Rates for Cabral with SGD Dils") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
R_Cabral_SGDDils
Plots for Respiration - Varari
- There is no significance at Varari for any environmental parameter or broad SGD dilutions
###################
## SGD dils
####################
### confirm a linear model
R_Varari_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="R") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(R_Varari_SGD_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(SGD_number + 1) 0.0066251 0.0066251 1 59.962 2.1962 0.1436### confirm a nonlinear model
R_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="R") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(R_Varari_SGD_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.012828 0.0064138 2 58.984 2.1691 0.1233R_Varari_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="R") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=SGD_number+1,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
# geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
# geom_hline(yintercept = 0) +
labs(x = "log SGD Dils",
y = "Rate (umol O2 cm2 hr-1)",
title = "R Rates for Varari with SGD Dils") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
R_Varari_SGDDils
Plots for Calcification - Cabral:
- Cabral was allegedly significant for temperature (0.01) and pH (0.01) with a nonlinear relationship - but not seeing that here for temp, only for pH
- Cabral is NOT significant for SGD Dils
###################
## pH
####################
### confirm nonlinear model
C_Cabral_pH_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(new_pH+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(C_Cabral_pH_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(new_pH + 1), 2) 0.043887 0.021943 2 47.288 3.2971 0.04566 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1C_Cabral_pH <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=new_pH,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
# geom_smooth(method="lm", formula = "y~log(x)", color="black") +
geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
# geom_hline(yintercept = 0) +
labs(x = "log pH",
y = "Rate (umol O2 cm2 hr-1)",
title = "C Rates for Cabral with pH") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
C_Cabral_pH
###################
## SGD dils
####################
### confirm a linear model
C_Cabral_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(C_Cabral_SGD_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(SGD_number + 1) 0.015433 0.015433 1 46.137 2.1383 0.1504### confirm a nonlinear model
C_Cabral_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(C_Cabral_SGD_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.01618 0.0080902 2 45.126 1.0993 0.3418C_Cabral_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Cabral") %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=SGD_number+1,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
#geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
# geom_hline(yintercept = 0) +
labs(x = "log SGD Dils",
y = "Rate (umol O2 cm2 hr-1)",
title = "C Rates for Cabral with SGD Dils") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
C_Cabral_SGDDils
Plots for Calcification - Varari
- Varari has a significant linear and nonlinear relationship with phosphate (<0.0001) and (0.0002) and TA (0.001)
- Varari has a significant … relationship with SGD
###################
## Phosphate
####################
### confirm nonlinear model
C_Varari_phosphate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(phosphate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="C") %>%
filter(phosphate_umolL<1) %>%
filter(!is.infinite(umol.cm2.hr)))
anova(C_Varari_phosphate_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## poly(log(phosphate_umolL + 1), 2) 0.011868 0.0059339 2 58.832 1.6788
## Pr(>F)
## poly(log(phosphate_umolL + 1), 2) 0.1954C_Varari_Phosphate <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="C") %>%
filter(phosphate_umolL<1) %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=phosphate_umolL,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
# geom_smooth(method="lm", formula = "y~log(x)", color="black") +
geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
# geom_hline(yintercept = 0) +
labs(x = "log phosphate",
y = "Rate (umol O2 cm2 hr-1)",
title = "C Rates for Varari with Phosphate") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
C_Varari_Phosphate
###################
## TA
####################
### confirm a nonlinear model
C_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(TA_initial+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(TA_initial<2550) %>% #### IF CHOOSE TO REMOVE HIGH POINTS, ONLY SIGNIFICANT IN NONLINEAR (otherwise significant in linear too)
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(C_Varari_SGD_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(TA_initial + 1), 2) 0.056876 0.028438 2 59.115 5.1362 0.008776 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1C_Varari_TA <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="C") %>%
# filter(TA_initial<2550) %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=TA_initial+1,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
geom_hline(yintercept = 0) +
labs(x = "log TA",
y = "Rate (umol O2 cm2 hr-1)",
title = "C Rates for Varari with TA (Minus the potential outliers)") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
C_Varari_TA
###################
## Silicate
####################
### nonlinear model
C_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(silicate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(C_Varari_SGD_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(silicate_umolL + 1), 2) 0.057443 0.028721 2 60.29 4.7675 0.01196
##
## poly(log(silicate_umolL + 1), 2) *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## not poly
C_Varari_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL+1)*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(C_Varari_silicate_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(silicate_umolL + 1) 0.046571 0.046571 1 107.262 6.9107 0.009826
## site 0.000212 0.000212 1 47.521 0.0315 0.859918
## log(silicate_umolL + 1):site 0.002363 0.002363 1 107.262 0.3507 0.554987
##
## log(silicate_umolL + 1) **
## site
## log(silicate_umolL + 1):site
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1C_Varari_silicate <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C") %>%
# filter(TA_initial<2550) %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=silicate_umolL+1,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=site)) +
scale_color_manual(values=pnw_palette("Bay", n=2),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~log(x)", color="black") +
geom_hline(yintercept = 0) +
labs(x = "log Silicate",
y = "Rate (umol O2 cm2 hr-1)",
title = "C Rates for Varari with Silicate") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
C_Varari_silicate
###################
## SGD dils
####################
### confirm a linear model
C_Varari_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(C_Varari_SGD_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(SGD_number + 1) 0.041992 0.041992 1 60.038 6.7823 0.01159 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### confirm a nonlinear model
C_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(C_Varari_SGD_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.082704 0.041352 2 59.049 7.3684 0.001391 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1C_Varari_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(site=="Varari") %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=SGD_number+1,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
scale_color_manual(values=pnw_palette("Bay", n=9),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~log(x)", color="black") +
geom_hline(yintercept = 0) +
labs(x = "log SGD Dils",
y = "Rate (umol O2 cm2 hr-1)",
title = "NEC Rates for Varari with SGD Dils") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
C_Varari_SGDDils
Plot nutrient levels with SGD dils to see how different the nutrients actually are across the dilutions
- just using as a check/baseline
Plot how many corals are net calcifying vs net dissolving
- this is supplementarl figure used as additional evidence to show difference in Varari gw vs Cabral
- Results: more fragments in Varari were net dissolving compared to Cabral (and at lower dilutions)
### should make a stacked bar plot with all corals in one group of dilution - % of which of those are positive or negative
###################################
## attempt 2
###################################
RespoR_Normalized_Full_withNEC_nutrients_calcifplot <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
select(!c("phosphate_umolL":"TA_initial", "salinity", "date")) %>%
mutate(rate_category = ifelse(umol.cm2.hr > 0, "positive", "negative")) %>%
group_by(SGD_number, site) %>%
dplyr::summarize(positive_count = sum(rate_category == "positive"),
negative_count = sum(rate_category == "negative"))
RespoR_Normalized_Full_withNEC_nutrients_calcifplot2 <- RespoR_Normalized_Full_withNEC_nutrients_calcifplot %>%
pivot_longer(cols = c(positive_count, negative_count),
names_to = "rate_category",
values_to = "count") %>%
mutate(total_count = sum(count),
proportion = count / total_count)
## plot
Net_CalcvsDissolution2 <- RespoR_Normalized_Full_withNEC_nutrients_calcifplot2 %>%
ggplot(aes(x=as.factor(SGD_number),
y=proportion,
fill=rate_category)) +
geom_bar(stat = "identity", position = "fill", width=0.4) +
#scale_fill_manual(values = ifelse (umol.cm2.hr > 0, "red", "black")) +
scale_fill_manual(values=c("red", "black"),
labels = c("positive_count" = "Net Calcifying", "negative_count" = "Net Dissolving")) +
theme_bw() +
# coord_trans(x="log") +
# geom_hline(yintercept = 0) +
facet_wrap(~site) +
#coord_trans(x ="log") +
labs(x = "% SGD by Volume",
y = "Proportion of Corals",
# title = "Net Calcifying vs Net Dissolving Corals per SGD Dilution",
fill = "Calcification State") +
theme(axis.text.x=element_text(size=14),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
Net_CalcvsDissolution2
ggsave(here::here("Outputs", "Results", "Final", "NetCalcifying_NetDissolving.jpg"),
width = 10, height=7)
########################
## statistics
#######################
CalcVsDissol <- lm(proportion ~ SGD_number*site, data = RespoR_Normalized_Full_withNEC_nutrients_calcifplot2)
anova(CalcVsDissol)## Analysis of Variance Table
##
## Response: proportion
## Df Sum Sq Mean Sq F value Pr(>F)
## SGD_number 1 0.00000 0.000000 0.0000 1.0000
## site 1 0.03840 0.038405 1.0651 0.3098
## SGD_number:site 1 0.00068 0.000684 0.0190 0.8913
## Residuals 32 1.15385 0.036058########################
## calculate percentage that are dissolving
#######################
percentDissolve <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C",
site=="Cabral") %>%
dplyr::mutate(percent_dissolved = (sum(umol.cm2.hr < 0) / n()) * 100)
##############################################
## look at raw data to see why dissolving is occurring at 0.01
##############################################
RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C",
SGD_number=="0.01") ### looks like it is only 3 colonies, 1 at Cabral (Col 5, but barely at -0.02), and 2 at Varari (Col 3 and 8)
RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C",
SGD_number=="0.03") ### only the same 2 colonies at Varari (Col 3 and 8)
RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C",
SGD_number=="0.05") ### only Varari Col 3
RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C",
SGD_number=="0.1") ### only Varari Col 3
RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C",
SGD_number=="0.5") ### Varari Col 3 and barely Col 8 (0.02)
RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C",
SGD_number=="1") ### Varari Col 1 and barely Col 3
RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C",
SGD_number=="2") ### one colony at Cabral (Col 10), Varari Col 3 and Col 8
RespoR_Normalized_Full_withNEC_nutrients_rawdata <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C",
SGD_number=="4") ### one colony at Cabral (Col 8), Varari Cols 2,3 and 4
###################################
## check GP rates of Varari Colony 3
###################################
RespoR_Normalized_Full_withNEC_nutrients_GPrawdata <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="GP",
new_colonynumber=="3"| new_colonynumber=="2")
RespoR_Normalized_Full_withNEC_nutrientsVGPonly <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="GP",
site=="Varari") %>%
filter(!is.infinite(umol.cm2.hr))
max(RespoR_Normalized_Full_withNEC_nutrientsVGPonly$umol.cm2.hr)## [1] 1.400944min(RespoR_Normalized_Full_withNEC_nutrientsVGPonly$umol.cm2.hr)## [1] 0.1810457###################################
## attempt 1
###################################
Net_CalcvsDissolution <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C") %>%
mutate(SGD_number = as.character(SGD_number),
rate_category = ifelse(umol.cm2.hr > 0, "positive", "negative")) %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=SGD_number,
y=umol.cm2.hr,
fill=rate_category)) +
geom_col(width = 0.8, colour = NA) +
#scale_fill_manual(values = ifelse (umol.cm2.hr > 0, "red", "black")) +
scale_fill_manual(values=c(positive='black',negative='red')) +
theme_bw() +
geom_hline(yintercept = 0) +
facet_wrap(~site) +
#coord_trans(x ="log") +
labs(x = "SGD values",
y = "Rate (umol O2 cm2 hr-1)",
title = "Net Calcifying vs Net Dissolving Corals per SGD Dilution") +
theme(axis.text.x=element_text(size=14),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
Net_CalcvsDissolution
Plot GP:R as a function of SGD
- Cabral is significant (0.025) but Varari is not
## First make new df
GP_RRatio_bySGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>%
dplyr::select(!"sample_ID") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
group_by(P_R, SGD_number, date, site, new_colonynumber) %>%
filter(P_R!="C", P_R!="NP") %>%
#summarize(GP_R = GP/R)
pivot_wider(names_from = P_R, values_from = umol.cm2.hr) %>%
mutate(GP_R = GP/R)
### plots
GP_R_bySGD_bothsites <- GP_RRatio_bySGDDils %>%
#filter(site=="Varari") %>%
mutate(new_colonynumber= as.factor(new_colonynumber)) %>%
ggplot(aes(x=(SGD_number+1),
y=GP_R)) +
geom_point(size=3, aes(color=new_colonynumber)) +
theme_bw() +
coord_trans(x="log") +
facet_wrap(~site, scales = "free_y") +
geom_smooth(method="lm", formula = "y~log(x)") +
labs(x = "SGD Dilutions",
y = "Gross Photosynthesis/Respiration Rates",
title = "SGD by GP:R Relationship for Varari") +
theme(axis.text.x=element_text(size=12),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
GP_R_bySGD_bothsites
## model testing
model_GP_R_SGD_Varari <- lmer(GP_R ~ log(SGD_number+1) + (1|new_colonynumber), data = GP_RRatio_bySGDDils %>%
filter(site=="Varari"))
anova(model_GP_R_SGD_Varari)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(SGD_number + 1) 0.26869 0.26869 1 60.003 3.6316 0.06148 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_GP_R_SGD_Cabral <- lmer(GP_R ~ log(SGD_number+1) + (1|new_colonynumber), data = GP_RRatio_bySGDDils %>%
filter(site=="Cabral"))
anova(model_GP_R_SGD_Cabral)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log(SGD_number + 1) 0.28731 0.28731 1 46.031 5.3534 0.0252 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NEW UPDATES FROM JUNE 19TH (MEETING WITH NYSSA):
Three main figures for publication will include:
Figure 1: Experimental Design, showing methods and set up
Figure 2: Correlation matrix (or PCA) by site to show the relationships of different parameters. Could be paired with a basic bar plot of end-member concentrations
Figure 3: 3 part figure including a plot for GP, R, and NEC with log(silicate) on the x axis and rates on the y axis, colored by site and only including geom_smooth line for NEC. GP can include a dashed line because Cabral was significant but Varari was not.
Exploratory figures: Net calcifying vs net dissolving, GP:C (potentially for index)
Figure 1: Conceptual (made in Biorender)
Figure 2: ReVamp’d (correlation matrix) (ignore these)
- this has all the significant correlations and labeled for each site
###################################
## Varari
###################################
## data
#corr_matrix_Varari <- RespoR_Normalized_Full_withNEC_nutrients %>%
# mutate(N_P = NN_umolL/phosphate_umolL) %>%
# filter(P_R=="R") %>%
# mutate("Salinity"=salinity,
# "TA" = TA_initial,
# "pH" = new_pH,
# "N+N" = NN_umolL,
# "Phosphate" = phosphate_umolL,
# "Silicate" = silicate_umolL) %>%
#"N+N:P" = N_P) %>%
# filter(site=="Varari") %>%
# select(!c("salinity", "silicate_umolL", "TA_initial", "new_pH", "N_P", "phosphate_umolL", "NN_umolL", #"umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL"))
# Compute a correlation matrix
#corr_Varari <- round(cor(corr_matrix_Varari,
# method="spearman"), 1)
# Compute a matrix of correlation p-values
#p.mat_Varari <- cor_pmat(corr_Varari)
# Visualize the correlation matrix
#ggcorrplot(corr_Varari, method = "circle")
# Reordering the correlation matrix using hierarchical clustering
#ggcorrplot(corr_Varari, hc.order = TRUE, outline.col = "white", method = "circle")
# Types of correlogram layout - Get the lower triangle
#ggcorrplot(corr_Varari, hc.order = TRUE, type = "lower", outline.col = "white", method = "circle")
# Change colors and theme
# Argument colors
#Varari_corr_matrix_plot <- ggcorrplot(corr_Varari,
# hc.order = TRUE, ## reordering the correlation matrix with hierarchial clustering
# type = "lower", ## correlogram layout - for the lower triangle
# outline.col = "white",
# # method = "circle", ## determines the shape
# ggtheme = ggplot2::theme_classic, ## changes the theme
# colors = c("steelblue3", "white", "tomato3"), ## determines color scheme
# sig.level = 0.05, ## default p-value
# lab = TRUE, ## adds correlation coefficient
# p.mat = p.mat_Varari, ## adds significance
# insig = "blank", ## makes insig values blank
# title = "Pearson Correlation Matrix for Environmental Parameters at Varari",
# legend.title = "Correlation",
# tl.srt = 45, ## slant of x axis
# tl.cex = 15, ## size of text axes
# lab_size = 5) +
#theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"),
# legend.title = element_text(size = 10, face="bold"),
# axis.text.x = element_text(size = 12, face="bold", color="black"),
# axis.text.y = element_text(size = 12, face="bold", color="black"))
#Varari_corr_matrix_plot
#ggsave(here::here("Outputs", "Results", "REDO_VarariCorrelation.jpg"),
# width=10, height=7)
###################################
## Cabral
###################################
## data
#corr_matrix_Cabral <- RespoR_Normalized_Full_withNEC_nutrients %>%
# mutate(N_P = NN_umolL/phosphate_umolL) %>%
# filter(P_R=="R") %>%
# mutate("Salinity"=salinity,
# "TA" = TA_initial,
# "pH" = new_pH,
# "N+N" = NN_umolL,
# "Phosphate" = phosphate_umolL,
# "Silicate" = silicate_umolL) %>%
# "N:P" = N_P) %>%
# filter(site=="Cabral") %>%
# select(!c("salinity", "silicate_umolL", "TA_initial", "new_pH", "N_P", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", # "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL"))
# Compute a correlation matrix
# corr_Cabral <- round(cor(corr_matrix_Cabral,
# method="spearman"), 1)
# Compute a matrix of correlation p-values
# p.mat_Cabral <- cor_pmat(corr_Cabral)
# Visualize the correlation matrix
# ggcorrplot(corr_Cabral, method = "circle")
#
# Change colors and theme
# Argument colors
# Cabral_corr_matrix_plot <- ggcorrplot(corr_Cabral,
# hc.order = TRUE, ## reordering the correlation matrix with hierarchial clustering
# type = "lower", ## correlogram layout - for the lower triangle
# outline.col = "white",
# method = "circle", ## determines the shape
# ggtheme = ggplot2::theme_classic, ## changes the theme
# colors = c("steelblue3", "white", "tomato3"), ## determines color scheme
#sig.level = 0.05, ## default p-value
# lab = TRUE, ## adds correlation coefficient
# p.mat = p.mat_Cabral, ## adds significance
# insig = "blank", ## makes insig values blank
# title = "Pearson Correlation Matrix for Environmental Parameters at Cabral",
# legend.title = "Correlation",
# tl.srt = 45, ## slant of x axis
# tl.cex = 15, ## size of text axes
# lab_size = 5) +
# theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"),
# # legend.title = element_text(size = 10, face="bold"),
# axis.text.x = element_text(size = 12, face="bold", color="black"),
# axis.text.y = element_text(size = 12, face="bold", color="black")) +
# labs(x="Location",
# y = expression(Phosphate~(mu*mol~L^-1)))
# Cabral_corr_matrix_plot
# ggsave(here::here("Outputs", "Results", "REDO_CabralCorrelation.jpg"),
# width=10, height=7)Figure 2: ReVamp’d again (ignore these)
- this is revamped but still includes all the squares further interpretation
###################################
## Varari
###################################
## data from above
# Compute a correlation matrix
# corr_Varari2 <- round(cor(corr_matrix_Varari), 1)
# Compute a matrix of correlation p-values
# p.mat_Varari2 <- cor_pmat(corr_Varari2)
# Specify the desired order of parameters
#param_orderV <- c("Salinity", "TA", "pH", "N+N", "Phosphate", "Silicate") # Replace with your actual parameter names
# Reorder the correlation matrix
#corr_Varari <- corr_Varari[param_orderV, param_orderV]
# Visualize the correlation matrix
# ggcorrplot(corr_Varari, method = "circle")
# Varari_corr_matrix_plot2 <- ggcorrplot(corr_Varari,
# hc.order = TRUE, ## reordering the correlation matrix with hierarchial clustering
# type = "lower", ## correlogram layout - for the lower triangle
# outline.col = "white",
# method = "circle", ## determines the shape
# ggtheme = ggplot2::theme_classic, ## changes the theme
# colors = c("steelblue3", "white", "tomato3"), ## determines color scheme
#sig.level = 0.05, ## default p-value
# lab = TRUE, ## adds correlation coefficient
#p.mat = p.mat_Varari, ## puts an X on the insignificant ones
# insig = "pch", ## makes insig values blank
#insig = label_sig,
# title = "Pearson Correlation Matrix for Environmental Parameters at Varari",
# legend.title = "Correlation",
# tl.srt = 45, ## slant of x axis
# # # tl.cex = 15, ## size of text axes
# lab_size = 5) +
#labs(y="Varari") +
# theme(legend.title = element_text(size = 10, face="bold"),
# axis.text.x = element_text(size = 12, face="bold", color="black"),
# axis.text.y = element_text(size = 12, face="bold", color="black"))
# Varari_corr_matrix_plot2
# ggsave(here::here("Outputs", "Results", "REDO_VarariCorrelation_extended.jpg"),
# width=10, height=7)
# #
###################################
## Cabral
###################################
# Compute a correlation matrix
# corr_Cabral2 <- round(cor(corr_matrix_Cabral), 1)
# Compute a matrix of correlation p-values
# p.mat_Cabral2 <- cor_pmat(corr_Cabral2)
# convert to long format so can add asterisk
#cor_long_Cabral <- melt(corr_Cabral2)
#p_long_Cabral <- melt(p.mat_Cabral2)
# Mark significant correlations with asterisks
#cor_long_Cabral$signif <- ifelse(p_long_Cabral$value < 0.05, "bold", "plain")
# Specify the desired order of parameters
# param_orderC <- c("Salinity", "TA", "pH", "N+N", "Phosphate", "Silicate") # Replace with your actual parameter names
# Reorder the correlation matrix
# corr_Cabral2 <- corr_Cabral2[param_orderC, param_orderC]
# ggcorrplot(corr_Cabral2, method = "circle")
# Change colors and theme
# Argument colors
# Cabral_corr_matrix_plot2 <- ggcorrplot(corr_Cabral2,
# hc.order = TRUE, ## reordering the correlation matrix with hierarchial clustering
# type = "lower", ## correlogram layout - for the lower triangle
# outline.col = "white",
# method = "circle", ## determines the shape
# ggtheme = ggplot2::theme_classic, ## changes the theme
# colors = c("steelblue3", "white", "tomato3"), ## determines color scheme
#sig.level = 0.05, ## default p-value
# lab = TRUE, ## adds correlation coefficient
# p.mat = p.mat_Cabral2, ## adds significance
#insig = "blank", ## makes insig values blank
# title = "Pearson Correlation Matrix for Environmental Parameters at Cabral",
# legend.title = "Correlation",
# tl.srt = 45, ## slant of x axis
# tl.cex = 15, ## size of text axes
# lab_size = 5) +
# theme(plot.title = element_text(hjust = 0.5, size=15, face="bold"),
# legend.title = element_text(size = 10, face="bold"),
# axis.text.x = element_text(size = 12, face="bold", color="black"),
# axis.text.y = element_text(size = 12, face="bold", color="black")) #+
#geom_text(data = cor_long_Cabral, aes(x = Var1, y = Var2, label = signif), color = "black", size = 5, inherit.aes = FALSE)
# Cabral_corr_matrix_plot2
# ggsave(here::here("Outputs", "Results", "REDO_CabralCorrelation_extended.jpg"),
# width=10, height=7)
###########################
## patch together
###########################
#### remove x-axis labels
# p_Ccorr <- Cabral_corr_matrix_plot2
# p_Vcorr <- Varari_corr_matrix_plot2
# combined_corr_plot <- p_Ccorr/p_Vcorr + plot_layout(guides = 'collect')
# combined_corr_plot +
# plot_annotation(tag_levels = 'A') &
# theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
# plot.tag = element_text(size = 18, face = "bold"))
# ggsave(here::here("Outputs", "Results", "Final", "CorrMatrices.jpg"),
# width=12, height=15)Figure 2a: correlation matrix ONLY with silicate AND spearman correlation - THIS IS THE ONE TO USE
- this is the final figure included for publication, which only has silicate compared to other env params AND uses the Spearman correlation which is better for nonlinear, monotonic data
###################################
## Varari
###################################
## data
corr_matrix_Varari_silicateonly <- RespoR_Normalized_Full_withNEC_nutrients %>%
# mutate(salinity = if_else(date == "6/17/23", salinity - 1, salinity), ### correcting for the weird day of salinity on 6/17 that was way higher than any other days. See code in Salinity_pH_TA_TestingRelationships for the deduction/evidence. Edited on 03/21/2025
mutate(N_P = NN_umolL/phosphate_umolL) %>%
filter(P_R=="C") %>%
mutate("Salinity"=salinity,
"TA" = TA_initial,
"pH" = new_pH,
"N+N" = NN_umolL,
"Phosphate" = phosphate_umolL,
"Silicate" = silicate_umolL) %>%
#"N+N:P" = N_P) %>%
filter(site=="Varari") %>%
select(!c("salinity", "silicate_umolL", "log_silicate", "log_silicate2", "TA_initial", "new_pH", "N_P", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL"))
## version 2.0
corr_Varari_spearman <- round(cor(corr_matrix_Varari_silicateonly, method="spearman"), 1)
# Compute a matrix of correlation p-values
p.mat_Varari_spearman <- cor_pmat(corr_Varari_spearman)
# Extract the correlations involving "Silicate"
corr_silicate_V <- corr_Varari_spearman[, "Silicate", drop = FALSE]
# Transpose the matrix to make "Silicate" the row instead of the column
corr_silicate_matrix_V <- t(as.matrix(corr_silicate_V))
# Convert the matrix to a long format
corr_melted_V <- melt(corr_silicate_matrix_V)
# Plot the correlation matrix vertically
Varari_corr_matrix_spearman <- ggplot(corr_melted_V, aes(x=Var2, y=Var1, fill=value)) +
geom_tile(color = "white") + # Adjust width of the tiles to make the plot narrower
scale_fill_gradient2(low = "steelblue3", high = "tomato3", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Spearman\nCorrelation") +
geom_text(aes(label = round(value, 2)), color = "black", size = 5) +
theme_minimal() +
theme(axis.text.x = element_text(vjust = 0.5, size = 12, face = "bold", color="black"),
axis.text.y = element_text(size = 12, face = "bold", color="black")) +
coord_flip() + # Flip coordinates for vertical layout
labs(x = "", y = "")
# Display the plot
Varari_corr_matrix_spearman
ggsave(here::here("Outputs", "Results", "Final", "Vsilicate_corrmatrix.jpg"),
width=7, height=12)
###################################
## Cabral
###################################
## data
corr_matrix_Cabral <- RespoR_Normalized_Full_withNEC_nutrients %>%
mutate(N_P = NN_umolL/phosphate_umolL) %>%
filter(P_R=="C") %>%
mutate("Salinity"=salinity,
"TA" = TA_initial,
"pH" = new_pH,
"N+N" = NN_umolL,
"Phosphate" = phosphate_umolL,
"Silicate" = silicate_umolL) %>%
# "N:P" = N_P) %>%
filter(site=="Cabral") %>%
select(!c("salinity", "log_silicate", "log_silicate2", "silicate_umolL", "TA_initial", "new_pH", "N_P", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL"))
# Compute a correlation matrix
corr_Cabral_spearman <- round(cor(corr_matrix_Cabral, method="spearman"), 1)
# Compute a matrix of correlation p-values
p.mat_Cabral_spearman <- cor_pmat(corr_Cabral_spearman)
# Specify the desired order of parameters
#param_orderC <- c("Salinity", "TA", "pH", "N+N", "Phosphate", "Silicate") # Replace with your actual parameter names
# Reorder the correlation matrix
#corr_Cabral2 <- corr_Cabral2[param_orderC, param_orderC]
# Extract the correlations involving "Silicate"
corr_silicate_C <- corr_Cabral_spearman[, "Silicate", drop = FALSE]
# Transpose the matrix to make "Silicate" the row instead of the column
corr_silicate_matrix_C <- t(as.matrix(corr_silicate_C))
# Convert the matrix to a long format
corr_melted_C <- melt(corr_silicate_matrix_C)
# Plot the correlation matrix vertically
Cabral_corr_matrix_spearman <- ggplot(corr_melted_C, aes(x=Var2, y=Var1, fill=value)) +
geom_tile(color = "white") + # Adjust width of the tiles to make the plot narrower
scale_fill_gradient2(low = "steelblue3", high = "tomato3", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Spearman\nCorrelation") +
geom_text(aes(label = round(value, 2)), color = "black", size = 5) +
theme_minimal() +
theme(axis.text.x = element_text(vjust = 0.5, size = 12, face = "bold", color="black"),
axis.text.y = element_text(size = 12, face = "bold", color="black")) +
coord_flip() + # Flip coordinates for vertical layout
labs(x = "", y = "")
Cabral_corr_matrix_spearman
ggsave(here::here("Outputs", "Results", "Final", "Csilicate_corrmatrix.jpg"),
width=7, height=12)
###########################
## patch together
###########################
#### remove x-axis labels
p_Ccorr_silicate <- Cabral_corr_matrix_spearman
p_Vcorr_silicate <- Varari_corr_matrix_spearman
combined_corr_plot_silicate <- p_Ccorr_silicate/p_Vcorr_silicate + plot_layout(guides = 'collect')
combined_corr_plot_silicate +
plot_annotation(tag_levels = 'A') &
theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
plot.tag = element_text(size = 18, face = "bold"))
ggsave(here::here("Outputs", "Results", "Final", "CorrMatrices.jpg"),
width=10, height=15)Getting p-values and r-squared coefficients for Spearman correlation plots
#################################
## Varari
#################################
# Calculate the Spearman correlation and p-values matrix
res_Varari_spearman <- rcorr(as.matrix(corr_matrix_Varari_silicateonly), type = "spearman")
# Extract the correlation matrix
corr_Varari_spearman <- res_Varari_spearman$r
corr_Varari_spearman## Salinity TA pH N+N Phosphate Silicate
## Salinity 1.00000000 0.06509682 -0.1481712 0.2119188 -0.3863537 -0.3078341
## TA 0.06509682 1.00000000 -0.5024418 0.7087955 0.3498518 0.8119126
## pH -0.14817116 -0.50244181 1.0000000 -0.5180828 -0.1335543 -0.4445489
## N+N 0.21191882 0.70879554 -0.5180828 1.0000000 0.2393305 0.5343213
## Phosphate -0.38635374 0.34985177 -0.1335543 0.2393305 1.0000000 0.4618148
## Silicate -0.30783408 0.81191260 -0.4445489 0.5343213 0.4618148 1.0000000# Extract the p-values matrix
p.mat_Varari_spearman <- res_Varari_spearman$P
p.mat_Varari_spearman## Salinity TA pH N+N Phosphate
## Salinity NA 5.896440e-01 2.175118e-01 7.603422e-02 8.749647e-04
## TA 0.5896439781 NA 8.015339e-06 4.647394e-12 2.783343e-03
## pH 0.2175117955 8.015339e-06 NA 3.702010e-06 2.668455e-01
## N+N 0.0760342177 4.647394e-12 3.702010e-06 NA 4.441475e-02
## Phosphate 0.0008749647 2.783343e-03 2.668455e-01 4.441475e-02 NA
## Silicate 0.0090129353 0.000000e+00 1.027994e-04 1.592406e-06 5.036852e-05
## Silicate
## Salinity 9.012935e-03
## TA 0.000000e+00
## pH 1.027994e-04
## N+N 1.592406e-06
## Phosphate 5.036852e-05
## Silicate NA# Calculate approximate r-squared by squaring the correlation coefficients
r_squared_Varari_spearman <- corr_Varari_spearman^2
r_squared_Varari_spearman## Salinity TA pH N+N Phosphate Silicate
## Salinity 1.000000000 0.004237596 0.02195469 0.04490959 0.14926921 0.09476182
## TA 0.004237596 1.000000000 0.25244777 0.50239112 0.12239626 0.65920206
## pH 0.021954692 0.252447770 1.00000000 0.26840979 0.01783675 0.19762370
## N+N 0.044909586 0.502391121 0.26840979 1.00000000 0.05727908 0.28549925
## Phosphate 0.149269214 0.122396261 0.01783675 0.05727908 1.00000000 0.21327287
## Silicate 0.094761820 0.659202062 0.19762370 0.28549925 0.21327287 1.00000000##########################
## only for silicate
# Extract correlations involving Silicate
corr_silicate_V <- corr_Varari_spearman[, "Silicate", drop = FALSE]
# Extract p-values involving Silicate
p_values_silicate_V <- p.mat_Varari_spearman[, "Silicate", drop = FALSE]
# Extract r-squared values involving Silicate
r_squared_silicate_V <- r_squared_Varari_spearman[, "Silicate", drop = FALSE]
# Print results for Silicate
print(corr_silicate_V)## Silicate
## Salinity -0.3078341
## TA 0.8119126
## pH -0.4445489
## N+N 0.5343213
## Phosphate 0.4618148
## Silicate 1.0000000print(p_values_silicate_V)## Silicate
## Salinity 9.012935e-03
## TA 0.000000e+00
## pH 1.027994e-04
## N+N 1.592406e-06
## Phosphate 5.036852e-05
## Silicate NAprint(r_squared_silicate_V)## Silicate
## Salinity 0.09476182
## TA 0.65920206
## pH 0.19762370
## N+N 0.28549925
## Phosphate 0.21327287
## Silicate 1.00000000# Identify significant correlations (p-value < 0.05)
significant_correlations <- p_values_silicate_V < 0.05
# Extract significant correlations
significant_corrs <- corr_silicate_V[significant_correlations]
# Extract the corresponding r-squared values
significant_r_squared <- r_squared_silicate_V[significant_correlations]
# Print significant correlations and their r-squared values
print(significant_corrs)## [1] -0.3078341 0.8119126 -0.4445489 0.5343213 0.4618148 NAprint(significant_r_squared)## [1] 0.09476182 0.65920206 0.19762370 0.28549925 0.21327287 NA#################################
## Cabral
#################################
# Calculate the Spearman correlation and p-values matrix
res_Cabral_spearman <- rcorr(as.matrix(corr_matrix_Cabral), type = "spearman")
# Extract the correlation matrix
corr_Cabral_spearman <- res_Cabral_spearman$r
corr_Cabral_spearman## Salinity TA pH N+N Phosphate
## Salinity 1.00000000 0.04362344 0.07363389 0.18414095 -0.5645513
## TA 0.04362344 1.00000000 0.02442002 -0.09770993 0.3998131
## pH 0.07363389 0.02442002 1.00000000 0.24030535 -0.1425907
## N+N 0.18414095 -0.09770993 0.24030535 1.00000000 -0.1536658
## Phosphate -0.56455133 0.39981309 -0.14259068 -0.15366585 1.0000000
## Silicate -0.80452665 0.03785682 0.14257365 -0.18964728 0.5856742
## Silicate
## Salinity -0.80452665
## TA 0.03785682
## pH 0.14257365
## N+N -0.18964728
## Phosphate 0.58567422
## Silicate 1.00000000# Extract the p-values matrix
p.mat_Cabral_spearman <- res_Cabral_spearman$P
p.mat_Cabral_spearman## Salinity TA pH N+N Phosphate
## Salinity NA 0.754118246 0.59669822 0.18255726 8.726575e-06
## TA 7.541182e-01 NA 0.86086202 0.48211665 2.741903e-03
## pH 5.966982e-01 0.860862023 NA 0.08006565 3.036804e-01
## N+N 1.825573e-01 0.482116652 0.08006565 NA 2.672614e-01
## Phosphate 8.726575e-06 0.002741903 0.30368044 0.26726141 NA
## Silicate 2.322587e-13 0.785792924 0.30373878 0.16959400 3.283811e-06
## Silicate
## Salinity 2.322587e-13
## TA 7.857929e-01
## pH 3.037388e-01
## N+N 1.695940e-01
## Phosphate 3.283811e-06
## Silicate NA# Calculate approximate r-squared by squaring the correlation coefficients
r_squared_Cabral_spearman <- corr_Cabral_spearman^2
r_squared_Cabral_spearman## Salinity TA pH N+N Phosphate
## Salinity 1.000000000 0.0019030045 0.0054219499 0.03390789 0.31871820
## TA 0.001903004 1.0000000000 0.0005963376 0.00954723 0.15985050
## pH 0.005421950 0.0005963376 1.0000000000 0.05774666 0.02033210
## N+N 0.033907888 0.0095472301 0.0577466635 1.00000000 0.02361319
## Phosphate 0.318718203 0.1598505036 0.0203321023 0.02361319 1.00000000
## Silicate 0.647263135 0.0014331385 0.0203272470 0.03596609 0.34301429
## Silicate
## Salinity 0.647263135
## TA 0.001433139
## pH 0.020327247
## N+N 0.035966090
## Phosphate 0.343014290
## Silicate 1.000000000##########################
## only for silicate
# Extract correlations involving Silicate
corr_silicate_C <- corr_Cabral_spearman[, "Silicate", drop = FALSE]
# Extract p-values involving Silicate
p_values_silicate_C <- p.mat_Cabral_spearman[, "Silicate", drop = FALSE]
# Extract r-squared values involving Silicate
r_squared_silicate_C <- r_squared_Cabral_spearman[, "Silicate", drop = FALSE]
# Print results for Silicate
print(corr_silicate_C)## Silicate
## Salinity -0.80452665
## TA 0.03785682
## pH 0.14257365
## N+N -0.18964728
## Phosphate 0.58567422
## Silicate 1.00000000print(p_values_silicate_C)## Silicate
## Salinity 2.322587e-13
## TA 7.857929e-01
## pH 3.037388e-01
## N+N 1.695940e-01
## Phosphate 3.283811e-06
## Silicate NAprint(r_squared_silicate_C)## Silicate
## Salinity 0.647263135
## TA 0.001433139
## pH 0.020327247
## N+N 0.035966090
## Phosphate 0.343014290
## Silicate 1.000000000# Identify significant correlations (p-value < 0.05)
significant_correlations <- p.mat_Cabral_spearman < 0.05
# Extract significant correlations
significant_corrs <- corr_Cabral_spearman[significant_correlations]
# Extract the corresponding r-squared values
significant_r_squared <- r_squared_Cabral_spearman[significant_correlations]
# Print significant correlations and their r-squared values
print(significant_corrs)## [1] NA -0.5645513 -0.8045267 NA 0.3998131 NA
## [7] NA -0.5645513 0.3998131 NA 0.5856742 -0.8045267
## [13] 0.5856742 NAprint(significant_r_squared)## [1] NA 0.3187182 0.6472631 NA 0.1598505 NA NA
## [8] 0.3187182 0.1598505 NA 0.3430143 0.6472631 0.3430143 NAFigure 2a (part2): Barplots for GW endmember data
- making bar graphs to illustrate differences between sites and from sites to offshore ambient SW
###########################################################
### making bar graphs to illustrate differences ###
##########################################################
#### SEE BACKGROUND SITE DATA SCRIPT FOR HOW OFFSHORE VALUES WERE CALCULATED ####
##############################
### plots for CSUNposium
###############################
sitecomparison_maxdata <- tribble(~Site, ~Parameter, ~Value,
"Varari", "pH", "7.73",
"Cabral", "pH", "6.89",
"Reef Crest", "pH", "7.9",
"Varari", "Nitrates + Nitrites (umol/L)", "42.49",
"Cabral", "Nitrates + Nitrites (umol/L)", "8.22",
"Reef Crest", "Nitrates + Nitrites (umol/L)", "0.58",
"Varari", "TA (umol/kg-1)", "5393.3",
"Cabral", "TA (umol/kg-1)", "1755.18",
"Reef Crest", "TA (umol/kg-1)", "2386.84",
"Varari", "Phosphate (umol/L)", "2.5",
"Cabral", "Phosphate (umol/L)", "3.61",
"Reef Crest", "Phosphate (umol/L)", "0.20",
"Varari", "Silicate (umol/L)", "725.61",
"Cabral", "Silicate (umol/L)", "855.01",
"Reef Crest", "Silicate (umol/L)", "3.68",
"Varari", "Salinity (psu)", "6.32",
"Cabral", "Salinity (psu)", "5.79",
"Reef Crest", "Salinity (psu)", "36.7") %>%
mutate(Site= as.factor(Site))
### for pH
SiteComparisons_pH <- sitecomparison_maxdata %>%
filter(Parameter=="pH") %>%
ggplot(aes(x=Site,
y=Value,
fill=Site)) +
geom_col() +
scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) +
scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"),
limits=c("Reef Crest", "Cabral", "Varari"),
guide="none") +
theme_classic() +
labs(x="Location",
y = expression(pH[T])) +
theme(axis.text.x=element_text(size=25),
axis.text.y=element_text(size=25),
axis.title.x=element_text(size=25),
axis.title.y=element_text(size=25))
SiteComparisons_pH
ggsave(here::here("Outputs", "Results", "SitevOffshore", "pH.jpg"))
## for TA
SiteComparisons_TA <- sitecomparison_maxdata %>%
filter(Parameter=="TA (umol/kg-1)") %>%
ggplot(aes(x=Site,
y=Value,
fill=Site)) +
geom_col() +
scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) +
scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"),
limits=c("Reef Crest", "Cabral", "Varari"),
guide="none") +
theme_classic() +
labs(x="Location",
y = expression(TA~(mu*mol~kg^-1))) +
theme(axis.text.x=element_text(size=25),
axis.text.y=element_text(size=25),
axis.title.x=element_text(size=25),
axis.title.y=element_text(size=25))
SiteComparisons_TA
ggsave(here::here("Outputs", "Results", "SitevOffshore", "TA.jpg"))
## Salinity
SiteComparisons_Salinity <- sitecomparison_maxdata %>%
filter(Parameter=="Salinity (psu)") %>%
ggplot(aes(x=Site,
y=Value,
fill=Site)) +
geom_col() +
scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) +
scale_y_discrete(limits = c("5.79", "6.32", "36.7")) +
scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"),
limits=c("Reef Crest", "Cabral", "Varari"),
guide="none") +
theme_classic() +
labs(x="Location",
y = "Salinity (psu)") +
theme(axis.text.x=element_text(size=25),
axis.text.y=element_text(size=25),
axis.title.x=element_text(size=25),
axis.title.y=element_text(size=25))
SiteComparisons_Salinity
ggsave(here::here("Outputs", "Results", "SitevOffshore", "Salinity.jpg"))
## for Phosphates
SiteComparisons_Phosphate <- sitecomparison_maxdata %>%
filter(Parameter=="Phosphate (umol/L)") %>%
ggplot(aes(x=Site,
y=Value,
fill=Site)) +
geom_col() +
scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) +
scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"),
limits=c("Reef Crest", "Cabral", "Varari"),
guide="none") +
theme_classic() +
labs(x="Location",
y = expression(Phosphate~(mu*mol~L^-1))) +
theme(axis.text.x=element_text(size=25),
axis.text.y=element_text(size=25),
axis.title.x=element_text(size=25),
axis.title.y=element_text(size=25))
SiteComparisons_Phosphate
ggsave(here::here("Outputs", "Results", "SitevOffshore", "Phosphate.jpg"))
## for Silicate
SiteComparisons_Silicate <- sitecomparison_maxdata %>%
filter(Parameter=="Silicate (umol/L)") %>%
ggplot(aes(x=Site,
y=Value,
fill=Site)) +
geom_col() +
scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) +
scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"),
limits=c("Reef Crest", "Cabral", "Varari"),
guide="none") +
theme_classic() +
labs(x="Location",
y = expression(Silicate~(mu*mol~L^-1))) +
theme(axis.text.x=element_text(size=25),
axis.text.y=element_text(size=25),
axis.title.x=element_text(size=25),
axis.title.y=element_text(size=25))
SiteComparisons_Silicate
ggsave(here::here("Outputs", "Results", "SitevOffshore", "Silicate.jpg"))
## nitrates
SiteComparisons_Nitrates <- sitecomparison_maxdata %>%
filter(Parameter=="Nitrates + Nitrites (umol/L)") %>%
ggplot(aes(x=Site,
y=as.numeric(Value),
fill=Site)) +
geom_col() +
scale_x_discrete(limits = c("Reef Crest", "Cabral", "Varari")) +
scale_y_continuous(trans='log10') +
#scale_y_continuous(trans = "log10", # Apply logarithmic transformation
# breaks = c(0.58, 8.22, 42.49), # Set specific breaks
# labels = c("0.58", "8.22", "42.49")) + # Custom labels for clarity
scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod2", "firebrick4"),
limits=c("Reef Crest", "Cabral", "Varari"),
guide="none") +
theme_classic() +
#scale_y_continuous(trans='log10') +
labs(y = expression(Nitrate+Nitrite~(mu*mol~L^-1))) +
theme(axis.text.x=element_text(size=25),
axis.text.y=element_text(size=25),
axis.title.x=element_text(size=25),
axis.title.y=element_text(size=25))
SiteComparisons_Nitrates
ggsave(here::here("Outputs", "Results", "SitevOffshore", "NN.jpg"))
### patch them all together
p1.1 <- SiteComparisons_Nitrates + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank())
p1.2 <- SiteComparisons_Phosphate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank())
p1.3 <- SiteComparisons_Silicate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank())
p1.4 <- SiteComparisons_Salinity + theme(axis.title.x = element_blank())
p1.5<- SiteComparisons_TA + theme(axis.title.x = element_blank())
p1.6 <- SiteComparisons_pH + theme(axis.title.x = element_blank())
combined_endmemberSW_plot <- (p1.1 + p1.2 + p1.3) / (p1.4 + p1.5 + p1.6) +
plot_layout(guides = 'collect') +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 24, face = "bold"))
combined_endmemberSW_plot
ggsave(here::here("Outputs", "Results", "Final", "rawGW_barplots.jpg"),
width=20, height = 12)
#plot.margin = unit(c(1, 1, 1, 1), "cm"), Figure 2b: Background information (plots with silicate vs parameters)
- for personal understanding - not to include in paper
##########################
## site data
##########################
background_site_plots <- RespoR_Normalized_Full_withNEC_nutrients %>%
# mutate(N_P = NN_umolL/phosphate_umolL) %>%
filter(P_R=="C") %>%
mutate("Salinity"=salinity,
"TA" = TA_initial,
"pH" = new_pH,
"NN" = NN_umolL,
"Phosphate" = phosphate_umolL,
"Silicate" = silicate_umolL) %>%
select(!c("salinity", "log_silicate", "log_silicate2", "silicate_umolL", "TA_initial", "new_pH", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "date", "P_R", "ammonia_umolL")) %>%
filter(Phosphate<0.9 & TA < 2550)
##########################
## ggplots - silciate vs everything else
##########################
### Phos
silicate_vs_phos <- background_site_plots %>%
ggplot(aes(x=Silicate,
y=Phosphate,
color=site)) +
geom_point(size=3) +
facet_wrap(~site, scales = "free_y") ## OUTLIER AT 0.9 for PHOS (removed)
silicate_vs_phos
### NN
silicate_vs_NN <- background_site_plots %>%
ggplot(aes(x=Silicate,
y= NN,
color=site)) +
geom_point(size=3) +
facet_wrap(~site, scales = "free_y")
silicate_vs_NN
### TA
silicate_vs_TA <- background_site_plots %>%
ggplot(aes(x= Silicate,
y= TA,
color=site)) +
geom_point(size=3) +
facet_wrap(~site, scales = "free_y") ## OUTLIER AT VARARI ABOVE 2550
silicate_vs_TA
### Salinity
silicate_vs_salinity <- background_site_plots %>%
ggplot(aes(x=Silicate,
y= Salinity,
color=site)) +
geom_point(size=3) +
facet_wrap(~site, scales = "free_y")
silicate_vs_salinity
### pH
silicate_vs_pH <- background_site_plots %>%
ggplot(aes(x=Silicate,
y= pH,
color=site)) +
geom_point(size=3) +
facet_wrap(~site, scales = "free_y")
silicate_vs_pH
############ patch together
silicate_by_otherparams <- silicate_vs_pH + silicate_vs_salinity + silicate_vs_TA + silicate_vs_NN + silicate_vs_phos
silicate_by_otherparams
ggsave(here::here("Outputs", "Results", "Final", "silicate_vs_otherparams.jpg"),
width=24, height = 12)Figure 3: Physio rate plots
### Respiration
R_silicate <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="R") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=silicate_umolL,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=site)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
theme_classic() +
coord_trans(x ="log") +
geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R == "R", !is.infinite(umol.cm2.hr)),
aes(x = silicate_umolL, y = umol.cm2.hr),
method = "lm", formula = y~poly(log(x),2), color = "black") +
labs(y = expression(Respiration~Rate~(mu*mol~O[2]~cm^-2~hr^-1)),
color = "Site") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=15, face="bold"),
axis.title.y=element_text(size=15, face="bold"),
legend.title = element_text(size=18,
face="bold"),
legend.text = element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20,
face="bold"))
R_silicate
ggsave(here::here("Outputs", "Results", "Final", "R_silicate.jpg"),
width=10, height=7)
RespoR_Normalized_Full_withNEC_nutrientsR <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C") %>%
filter(site=="Cabral") %>%
filter(!is.infinite(umol.cm2.hr))
max(RespoR_Normalized_Full_withNEC_nutrientsR$umol.cm2.hr)## [1] 0.4004983min(RespoR_Normalized_Full_withNEC_nutrientsR$umol.cm2.hr)## [1] -0.06829634###############################
## Gross Photosynthesis
###############################
GP_silicate <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=silicate_umolL,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=site)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
theme_classic() +
coord_trans(x ="log") +
geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R == "GP", !is.infinite(umol.cm2.hr)),
aes(x = silicate_umolL, y = umol.cm2.hr),
method = "lm", formula = y~poly(log(x),2), color = "black") +
labs(x = expression(Silicate~(mu*mol~L^-1)),
y = expression(Gross~Photosynthesis~Rate~(mu*mol~O[2]~cm^-2~hr^-1)),
color = "Site") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=20, face="bold"),
axis.title.y=element_text(size=20, face="bold"),
legend.title = element_text(size=18,
face="bold"),
legend.text = element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20,
face="bold"))
GP_silicate
ggsave(here::here("Outputs", "Results", "GP_silicate.jpg"),
width=10, height=7)
###############################
## Calcification
###############################
C_silicate_ECRS_firebrick4 <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=silicate_umolL,
y=umol.cm2.hr)) +
geom_point(size=5, aes(color=site)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
theme_classic() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
#scale_x_discrete(breaks = c(1, 10, 20)) +
geom_hline(yintercept = 0) +
geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R == "C", !is.infinite(umol.cm2.hr)),
aes(x = silicate_umolL, y = umol.cm2.hr),
method = "lm", formula = y~log(x), color = "black") +
labs(x = expression(Silicate~(mu*mol~L^-1)),
y = expression(Calcification~Rate~(mu*mol~CaCO[3]~cm^-2~hr^-1)),
color = "Site") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=15, face="bold"),
axis.title.y=element_text(size=15, face="bold"),
legend.title = element_text(size=18,
face="bold"),
legend.text = element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20,
face="bold"))
C_silicate_ECRS_firebrick4
ggsave(here::here("Outputs", "Results", "C_silicate.jpg"),
width=10, height=7)
#############################################
### make these plots figure worthy using patchwork
#############################################
#### remove x-axis labels
p1 <- R_silicate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(),
legend.title = element_blank(), legend.text = element_blank(), legend.position = "none",
axis.title.y=element_text(size=15, face="bold")) +
scale_x_continuous(breaks = c(2, 10, 20))
p2 <- GP_silicate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(),
axis.title.y=element_text(size=15, face="bold")) +
scale_x_continuous(breaks = c(2, 10, 20))
p3 <- C_silicate_ECRS_firebrick4 +
theme(legend.title = element_blank(), legend.text = element_blank(), legend.position = "none",
axis.title.x=element_text(size=15, face="bold"),
axis.title.y=element_text(size=15, face="bold")) +
scale_x_continuous(breaks = c(2, 10, 20))
combined_column <- p1/p2/p3 + plot_layout(guides = 'collect')
combined_column +
plot_annotation(tag_levels = 'A') &
theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
plot.tag = element_text(size = 24, face = "bold"))
ggsave(here::here("Outputs", "Results", "Final", "CombinedRates.jpg"),
width=10, height=16)Figure 3: Re’Vampd with GP and R nonlinear
###################
### GP 2.0
###################
#### confirm models are significant for both sites
### confirm nonlinear model
GP_Varari_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(silicate_umolL+1),2)*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
#filter(site=="Varari") %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_Varari_silicate_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## poly(log(silicate_umolL + 1), 2) 0.210557 0.105279 2 106.556 6.2548
## site 0.009422 0.009422 1 11.731 0.5598
## poly(log(silicate_umolL + 1), 2):site 0.045241 0.022621 2 106.556 1.3439
## Pr(>F)
## poly(log(silicate_umolL + 1), 2) 0.002701 **
## site 0.469077
## poly(log(silicate_umolL + 1), 2):site 0.265202
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(GP_Varari_silicate_nonlinearmodel)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ poly(log(silicate_umolL + 1), 2) * site + (1 |
## new_colonynumber)
## Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "GP") %>%
## filter(!is.infinite(umol.cm2.hr))
##
## REML criterion at convergence: -120.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9126 -0.6755 0.1116 0.6361 3.4111
##
## Random effects:
## Groups Name Variance Std.Dev.
## new_colonynumber (Intercept) 0.00760 0.08718
## Residual 0.01683 0.12974
## Number of obs: 122, groups: new_colonynumber, 14
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 0.76218 0.03994 11.72728
## poly(log(silicate_umolL + 1), 2)1 -0.53415 0.22012 104.47568
## poly(log(silicate_umolL + 1), 2)2 0.43175 0.23982 105.71662
## siteVarari -0.03954 0.05285 11.73145
## poly(log(silicate_umolL + 1), 2)1:siteVarari 0.43517 0.28499 105.40538
## poly(log(silicate_umolL + 1), 2)2:siteVarari 0.08519 0.30228 108.77719
## t value Pr(>|t|)
## (Intercept) 19.081 3.42e-10 ***
## poly(log(silicate_umolL + 1), 2)1 -2.427 0.0169 *
## poly(log(silicate_umolL + 1), 2)2 1.800 0.0747 .
## siteVarari -0.748 0.4691
## poly(log(silicate_umolL + 1), 2)1:siteVarari 1.527 0.1298
## poly(log(silicate_umolL + 1), 2)2:siteVarari 0.282 0.7786
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl((_L+1),2)1 pl((_L+1),2)2 sitVrr p((_L+1),2)1:
## pl((_L+1),2)1 -0.068
## pl((_L+1),2)2 0.070 -0.413
## siteVarari -0.756 0.052 -0.053
## p((_L+1),2)1: 0.053 -0.772 0.319 -0.026
## p((_L+1),2)2: -0.056 0.328 -0.793 0.029 -0.199r.squaredGLMM(GP_Varari_silicate_nonlinearmodel)## R2m R2c
## [1,] 0.09963739 0.3797085#### now plot
GP_silicate2 <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=silicate_umolL,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=site)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
theme_classic() +
coord_trans(x ="log") +
geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R == "GP", site == "Varari", !is.infinite(umol.cm2.hr)),
aes(x = silicate_umolL, y = umol.cm2.hr),
method = "lm", formula = y~poly(log(x),2), color = "firebrick4") +
geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R == "GP", site == "Cabral", !is.infinite(umol.cm2.hr)),
aes(x = silicate_umolL, y = umol.cm2.hr),
method = "lm", formula = y~poly(log(x),2), color = "darkgoldenrod2") +
labs(x = expression(Silicate~(mu*mol~L^-1)),
y = expression(Gross~Photosynthesis~Rate~(mu*mol~O[2]~cm^2~hr^-1)),
color = "Site") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=20, face="bold"),
axis.title.y=element_text(size=20, face="bold"),
legend.title = element_text(size=18,
face="bold"),
legend.text = element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20,
face="bold"))
GP_silicate2 
###########################
## patch together
###########################
#### remove x-axis labels
p1 <- R_silicate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(),
legend.title = element_blank(), legend.text = element_blank(), legend.position = "none",
axis.title.y=element_text(size=15, face="bold")) +
scale_x_continuous(breaks = c(2, 10, 20))
p2 <- GP_silicate2 + theme(axis.title.x = element_blank(), axis.text.x = element_blank(),
axis.title.y=element_text(size=15, face="bold")) +
scale_x_continuous(breaks = c(2, 10, 20))
p3 <- C_silicate_ECRS_firebrick4 +
theme(legend.title = element_blank(), legend.text = element_blank(), legend.position = "none",
axis.title.x=element_text(size=15, face="bold"),
axis.title.y=element_text(size=15, face="bold")) +
scale_x_continuous(breaks = c(2, 10, 20))
combined_column <- p1/p2/p3 + plot_layout(guides = 'collect')
combined_column +
plot_annotation(tag_levels = 'A') &
theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
plot.tag = element_text(size = 18, face = "bold"))
ggsave(here::here("Outputs", "Results", "Final", "CombinedRates.jpg"),
width=10, height=15)FINAL MODELS USED IN RESULTS:
RESULTS SHOW THAT: - linear model is better fit for R, GP, and C - decided on ANCOVA design rather than analyzing sites separately – acknowledging that they are different BUT want to test for the relationship between the two to see HOW they are different and how the GP, R, and C rates are different in the fragments exposed to gw from both sites - model set-up: “model <- lmer(umol.cm2.hr ~ log_silicate*site + log_silicate2 + (1|new_colonynumber), data=data” – and then filter for specifically R, GP, or C
- Respiration: log_silicate, 2 (poly) = 0.0065 and log_silicate is 0.0045, so both significant. The interaction between silicate and site = 0.95 and just site (0.2) are both insignificant
- Photosynthesis: Poly (0.001) and linear silicate (0.0005) both significant. Interaction and just site are both insignificant
- Calcification: ONLY testing linear. Logged silicate = 0.01 so significant. Interaction and just site are both insignificant
RespoR_Normalized_Full_withNEC_nutrients <- RespoR_Normalized_Full_withNEC_nutrients %>%
mutate(log_silicate = log(silicate_umolL+1),
log_silicate2 = (log_silicate^2),
date=factor(date))
########################
## R
########################
### nonlinear model
#R_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + log_silicate2 + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
# filter(P_R=="R") %>%
#filter(!is.infinite(umol.cm2.hr)))
R_silicate_nonlinearmodel2 <- lmer(umol.cm2.hr ~ log_silicate*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(R_silicate_nonlinearmodel2)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_silicate 0.046571 0.046571 1 107.262 6.9107 0.009826 **
## site 0.000212 0.000212 1 47.521 0.0315 0.859918
## log_silicate:site 0.002363 0.002363 1 107.262 0.3507 0.554987
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(R_silicate_nonlinearmodel2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log_silicate * site + (1 | new_colonynumber)
## Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "C") %>%
## filter(!is.infinite(umol.cm2.hr))
##
## REML criterion at convergence: -214.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6824 -0.5040 -0.0080 0.5366 2.8061
##
## Random effects:
## Groups Name Variance Std.Dev.
## new_colonynumber (Intercept) 0.006216 0.07884
## Residual 0.006739 0.08209
## Number of obs: 122, groups: new_colonynumber, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.16360 0.05017 48.53078 3.261 0.00204 **
## log_silicate -0.02525 0.01844 106.08033 -1.369 0.17385
## siteVarari -0.01173 0.06611 47.52070 -0.177 0.85992
## log_silicate:siteVarari -0.01468 0.02480 107.26152 -0.592 0.55499
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_slc sitVrr
## log_silicat -0.733
## siteVarari -0.759 0.557
## lg_slct:stV 0.545 -0.744 -0.730check_model(R_silicate_nonlinearmodel2)
AIC(R_silicate_nonlinearmodel2)## [1] -202.3236r.squaredGLMM(R_silicate_nonlinearmodel2)## R2m R2c
## [1,] 0.05804063 0.5099971### linear model
#R_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL+1)*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
#filter(site=="Varari") %>%
# filter(P_R=="R") %>%
# filter(!is.infinite(umol.cm2.hr)))
R_bothsites_silicate_linearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="R") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(R_bothsites_silicate_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_silicate 0.0042170 0.0042170 1 107.372 1.2479 0.2664
## site 0.0018190 0.0018190 1 53.797 0.5383 0.4663
## log_silicate:site 0.0014839 0.0014839 1 107.372 0.4391 0.5090summary(R_bothsites_silicate_linearmodel)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log_silicate * site + (1 | new_colonynumber)
## Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "R") %>%
## filter(!is.infinite(umol.cm2.hr))
##
## REML criterion at convergence: -297.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.77505 -0.58754 0.08962 0.56703 2.36504
##
## Random effects:
## Groups Name Variance Std.Dev.
## new_colonynumber (Intercept) 0.002584 0.05083
## Residual 0.003379 0.05813
## Number of obs: 122, groups: new_colonynumber, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.361900 0.034255 55.124914 10.565 7.3e-15 ***
## log_silicate -0.003988 0.013059 105.976258 -0.305 0.761
## siteVarari -0.033095 0.045107 53.796820 -0.734 0.466
## log_silicate:siteVarari -0.011629 0.017549 107.372117 -0.663 0.509
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_slc sitVrr
## log_silicat -0.761
## siteVarari -0.759 0.578
## lg_slct:stV 0.566 -0.744 -0.757AIC(R_bothsites_silicate_linearmodel)## [1] -285.7517########################
## GP
########################
### nonlinear model
GP_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + log_silicate2 + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(GP_silicate_nonlinearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_silicate 0.212785 0.212785 1 112.396 12.7387 0.000528 ***
## site 0.051338 0.051338 1 80.138 3.0734 0.083404 .
## log_silicate2 0.185984 0.185984 1 111.720 11.1343 0.001151 **
## log_silicate:site 0.043913 0.043913 1 105.552 2.6289 0.107916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(GP_silicate_nonlinearmodel)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## umol.cm2.hr ~ log_silicate * site + log_silicate2 + (1 | new_colonynumber)
## Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "GP") %>%
## filter(!is.infinite(umol.cm2.hr))
##
## REML criterion at convergence: -110.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9533 -0.6610 0.1187 0.6322 3.4536
##
## Random effects:
## Groups Name Variance Std.Dev.
## new_colonynumber (Intercept) 0.007516 0.0867
## Residual 0.016704 0.1292
## Number of obs: 122, groups: new_colonynumber, 14
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.39873 0.17566 116.51375 7.963 1.24e-12 ***
## log_silicate -0.57275 0.15819 111.57997 -3.621 0.000444 ***
## siteVarari -0.16698 0.09525 80.13835 -1.753 0.083404 .
## log_silicate2 0.11436 0.03427 111.71969 3.337 0.001151 **
## log_silicate:siteVarari 0.06559 0.04045 105.55218 1.621 0.107916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_slc sitVrr lg_sl2
## log_silicat -0.962
## siteVarari -0.524 0.358
## log_silict2 0.917 -0.983 -0.250
## lg_slct:stV 0.485 -0.398 -0.834 0.271AIC(GP_silicate_nonlinearmodel)## [1] -96.5832r.squaredGLMM(GP_silicate_nonlinearmodel)## R2m R2c
## [1,] 0.09906217 0.3786475### linear model
#GP_silicate_linearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
#filter(site=="Varari") %>%
# filter(P_R=="GP") %>%
#filter(!is.infinite(umol.cm2.hr)))
#anova(GP_silicate_linearmodel)
#summary(GP_silicate_linearmodel)
#AIC(GP_silicate_linearmodel)
########################
## Calcification
########################
### nonlinear model
#C_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + log_silicate2 + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
#filter(site=="Varari") %>%
# filter(P_R=="C") %>%
# filter(!is.infinite(umol.cm2.hr)))
#anova(C_Varari_silicate_nonlinearmodel)
#summary(C_Varari_silicate_nonlinearmodel)
#AIC(C_Varari_silicate_nonlinearmodel)
### linear model
C_silicate_linearmodel <- lmer(umol.cm2.hr ~ log_silicate*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(C_silicate_linearmodel)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_silicate 0.046571 0.046571 1 107.262 6.9107 0.009826 **
## site 0.000212 0.000212 1 47.521 0.0315 0.859918
## log_silicate:site 0.002363 0.002363 1 107.262 0.3507 0.554987
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(C_Varari_silicate_linearmodel)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log(silicate_umolL + 1) * site + (1 | new_colonynumber)
## Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "C") %>%
## filter(!is.infinite(umol.cm2.hr))
##
## REML criterion at convergence: -214.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6824 -0.5040 -0.0080 0.5366 2.8061
##
## Random effects:
## Groups Name Variance Std.Dev.
## new_colonynumber (Intercept) 0.006216 0.07884
## Residual 0.006739 0.08209
## Number of obs: 122, groups: new_colonynumber, 14
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.16360 0.05017 48.53078 3.261
## log(silicate_umolL + 1) -0.02525 0.01844 106.08033 -1.369
## siteVarari -0.01173 0.06611 47.52070 -0.177
## log(silicate_umolL + 1):siteVarari -0.01468 0.02480 107.26152 -0.592
## Pr(>|t|)
## (Intercept) 0.00204 **
## log(silicate_umolL + 1) 0.17385
## siteVarari 0.85992
## log(silicate_umolL + 1):siteVarari 0.55499
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg(_L+1) sitVrr
## lg(slc_L+1) -0.733
## siteVarari -0.759 0.557
## lg(s_L+1):V 0.545 -0.744 -0.730AIC(C_Varari_silicate_linearmodel)## [1] -202.3236r.squaredGLMM(C_silicate_linearmodel)## R2m R2c
## [1,] 0.05804063 0.5099971########################################################
########################################################Models worked on August 29th, 2024
#####################
## Respiration
#####################
modeltest_R <- lmer(umol.cm2.hr ~ log_silicate*site + I(log_silicate^2) + (1|date/new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="R") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(modeltest_R)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_silicate 0.0290541 0.0290541 1 106.795 9.2865 0.002909 **
## site 0.0020049 0.0020049 1 9.596 0.6408 0.442786
## I(log_silicate^2) 0.0265717 0.0265717 1 106.625 8.4931 0.004345 **
## log_silicate:site 0.0000016 0.0000016 1 105.325 0.0005 0.981786
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(modeltest_R)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log_silicate * site + I(log_silicate^2) + (1 |
## date/new_colonynumber)
## Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "R") %>%
## filter(!is.infinite(umol.cm2.hr))
##
## REML criterion at convergence: -307.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.49031 -0.60178 0.00364 0.73007 2.41943
##
## Random effects:
## Groups Name Variance Std.Dev.
## new_colonynumber:date (Intercept) 0.0001486 0.01219
## date (Intercept) 0.0044872 0.06699
## Residual 0.0031286 0.05593
## Number of obs: 122, groups: new_colonynumber:date, 14; date, 7
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.683e-01 8.496e-02 6.160e+01 6.690 7.59e-09 ***
## log_silicate -2.032e-01 6.949e-02 1.066e+02 -2.924 0.00422 **
## siteVarari -5.035e-02 6.290e-02 9.596e+00 -0.801 0.44279
## I(log_silicate^2) 4.389e-02 1.506e-02 1.066e+02 2.914 0.00434 **
## log_silicate:siteVarari 4.011e-04 1.753e-02 1.053e+02 0.023 0.98179
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_slc sitVrr I(_^2)
## log_silicat -0.873
## siteVarari -0.548 0.231
## I(lg_slc^2) 0.833 -0.984 -0.162
## lg_slct:stV 0.428 -0.385 -0.545 0.260#####################
## Photosynthesis
#####################
modeltest_GP <- lmer(umol.cm2.hr ~ log_silicate*site + I(log_silicate^2) + (1|date/new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="GP") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(modeltest_GP)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_silicate 0.220744 0.220744 1 109.495 13.2207 0.0004236 ***
## site 0.050153 0.050153 1 36.857 3.0037 0.0914350 .
## I(log_silicate^2) 0.193149 0.193149 1 109.259 11.5680 0.0009381 ***
## log_silicate:site 0.042678 0.042678 1 105.478 2.5560 0.1128650
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(modeltest_GP)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log_silicate * site + I(log_silicate^2) + (1 |
## date/new_colonynumber)
## Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "GP") %>%
## filter(!is.infinite(umol.cm2.hr))
##
## REML criterion at convergence: -110.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9118 -0.6580 0.1101 0.6406 3.4806
##
## Random effects:
## Groups Name Variance Std.Dev.
## new_colonynumber:date (Intercept) 0.005917 0.07692
## date (Intercept) 0.001941 0.04406
## Residual 0.016697 0.12922
## Number of obs: 122, groups: new_colonynumber:date, 14; date, 7
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.41148 0.17736 107.13408 7.958 1.94e-12 ***
## log_silicate -0.58502 0.15884 109.21295 -3.683 0.000360 ***
## siteVarari -0.17129 0.09883 36.85738 -1.733 0.091435 .
## I(log_silicate^2) 0.11705 0.03442 109.25888 3.401 0.000938 ***
## log_silicate:siteVarari 0.06469 0.04046 105.47770 1.599 0.112865
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_slc sitVrr I(_^2)
## log_silicat -0.956
## siteVarari -0.523 0.345
## I(lg_slc^2) 0.912 -0.983 -0.242
## lg_slct:stV 0.478 -0.394 -0.803 0.268PLOTS BASED ON MODELS: RESPIRATION
###################
### R
###################
modeltest_R <- lmer(umol.cm2.hr ~ log_silicate*site + I(log_silicate^2) + (1|date/new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="R") %>%
filter(!is.infinite(umol.cm2.hr)))
anova(modeltest_R)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## log_silicate 0.0290541 0.0290541 1 106.795 9.2865 0.002909 **
## site 0.0020049 0.0020049 1 9.596 0.6408 0.442786
## I(log_silicate^2) 0.0265717 0.0265717 1 106.625 8.4931 0.004345 **
## log_silicate:site 0.0000016 0.0000016 1 105.325 0.0005 0.981786
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(modeltest_R)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: umol.cm2.hr ~ log_silicate * site + I(log_silicate^2) + (1 |
## date/new_colonynumber)
## Data: RespoR_Normalized_Full_withNEC_nutrients %>% filter(P_R == "R") %>%
## filter(!is.infinite(umol.cm2.hr))
##
## REML criterion at convergence: -307.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.49031 -0.60178 0.00364 0.73007 2.41943
##
## Random effects:
## Groups Name Variance Std.Dev.
## new_colonynumber:date (Intercept) 0.0001486 0.01219
## date (Intercept) 0.0044872 0.06699
## Residual 0.0031286 0.05593
## Number of obs: 122, groups: new_colonynumber:date, 14; date, 7
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.683e-01 8.496e-02 6.160e+01 6.690 7.59e-09 ***
## log_silicate -2.032e-01 6.949e-02 1.066e+02 -2.924 0.00422 **
## siteVarari -5.035e-02 6.290e-02 9.596e+00 -0.801 0.44279
## I(log_silicate^2) 4.389e-02 1.506e-02 1.066e+02 2.914 0.00434 **
## log_silicate:siteVarari 4.011e-04 1.753e-02 1.053e+02 0.023 0.98179
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) lg_slc sitVrr I(_^2)
## log_silicat -0.873
## siteVarari -0.548 0.231
## I(lg_slc^2) 0.833 -0.984 -0.162
## lg_slct:stV 0.428 -0.385 -0.545 0.260### based on reference from stack overflow on how to visualize a mixed model with random effects
### create a new dataframe that has x var (silicate), z var (color by colony number), and then create y values using predict function
### logged or not logged silicate?
### first isolate just R values
R_only <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(!is.infinite(umol.cm2.hr)) %>%
filter(P_R=="R")
write_csv(R_only, here::here("Data", "R_only.csv"))
### now make new df
plot_df_R <- data.frame(silicate = R_only$silicate_umolL,
pred_umol = predict(modeltest_R, newdata = R_only),
new_colonynumber = as.factor(R_only$new_colonynumber))
### define the desired color groups based on colony number (z variable)
plot_df_R$color_group <- ifelse(plot_df_R$new_colonynumber %in% as.character(1:8), "red2",
ifelse(plot_df_R$new_colonynumber %in% as.character(13:18), "gold2", NA))
### define the gradients (Varari in red, Cabral in blue)
Varari_gradient <- scales::gradient_n_pal(c("red2", "coral4"))(seq(0, 1, length.out = length(unique(plot_df_R$new_colonynumber[plot_df_R$color_group == "red2"]))))
Cabral_gradient <- scales::gradient_n_pal(c("gold2", "goldenrod"))(seq(0, 1, length.out = length(unique(plot_df_R$new_colonynumber[plot_df_R$color_group == "gold2"]))))
### Create the ggplot with conditional gradients
Rplot <- plot_df_R %>%
ggplot(aes(x = silicate,
y = pred_umol,
color = new_colonynumber)) +
geom_point(size = 3, alpha = 0.5,
data = R_only,
aes(x= silicate_umolL,
y= umol.cm2.hr,
color = as.factor(new_colonynumber))) +
geom_line(data = plot_df_R, linetype = 1) +
coord_trans(x ="log") +
# geom_abline(aes(intercept = fixef(modeltest_R)[1], slope = fixef(modeltest_R)[2],
# linetype = "Fixed effect"), key_glyph = draw_key_path,
# linewidth = 1) +
scale_color_manual(values = c(Varari_gradient, Cabral_gradient)) +
labs(linetype = NULL, color = "Colony Number", color_group = "Site") +
# annotation_custom(grid::textGrob(
#bquote(bold(~y ==
# .(round(fixef(modeltest_R)[1], 3)) + .(round(fixef(modeltest_R)[2], 3)) * x)),
# x = 0.05, y = 0.9, hjust = 0, gp = grid::gpar(cex = 1.2))) +
ggtitle("Illustration of mixed effects model with random intercept for RESPIRATION") +
theme(legend.position = "top",
plot.title = element_text(face = 2, hjust = 0.5),
panel.border = element_rect(fill = NA, linewidth = 0.2))
Rplot
PLOTS BASED ON MODELS: PHOTOSYNTHESIS
###################
### GP
###################
### based on reference from stack overflow on how to visualize a mixed model with random effects
### create a new dataframe that has x var (silicate), z var (color by colony number), and then create y values using predict function
### logged or not logged silicate?
### first isolate just R values
GP_only <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(!is.infinite(umol.cm2.hr)) %>%
filter(P_R=="GP")
### now make new df
plot_df_GP <- data.frame(silicate = GP_only$log_silicate,
pred_umol = predict(modeltest_GP, newdata = GP_only),
new_colonynumber = as.factor(GP_only$new_colonynumber))
### define the desired color groups based on colony number (z variable)
plot_df_GP$color_group <- ifelse(plot_df_GP$new_colonynumber %in% as.character(1:8), "firebrick4",
ifelse(plot_df_GP$new_colonynumber %in% as.character(13:18), "goldenrod", NA))
### define the gradients (Varari in red, Cabral in blue)
red_gradient <- scales::gradient_n_pal(c("firebrick4"))(seq(0, 1, length.out = length(unique(plot_df_GP$new_colonynumber[plot_df_GP$color_group == "firebrick4"]))))
blue_gradient <- scales::gradient_n_pal(c("goldenrod"))(seq(0, 1, length.out = length(unique(plot_df_GP$new_colonynumber[plot_df_GP$color_group == "goldenrod"]))))
### Create the ggplot with conditional gradients
GPplot <- plot_df_GP %>%
ggplot(aes(x = silicate,
y = pred_umol,
color = new_colonynumber)) +
geom_point(size = 3, alpha = 0.5,
data = GP_only,
aes(x= log_silicate,
y= umol.cm2.hr,
color = as.factor(new_colonynumber))) +
geom_line(data = plot_df_GP, linetype = 2) +
# geom_abline(aes(intercept = fixef(modeltest_GP)[1], slope = fixef(modeltest_GP)[2],
# linetype = "Fixed effect"), key_glyph = draw_key_path,
# linewidth = 1) +
scale_color_manual(values = c(red_gradient, blue_gradient)) +
labs(linetype = NULL, color = "Colony Number") +
annotation_custom(grid::textGrob(
bquote(bold(~y ==
.(round(fixef(modeltest_GP)[1], 3)) + .(round(fixef(modeltest_GP)[2], 3)) * x)),
x = 0.05, y = 0.9, hjust = 0, gp = grid::gpar(cex = 1.2))) +
ggtitle("Illustration of mixed effects model with random intercept for PHOTOSYNTHESIS") +
theme(legend.position = "top",
plot.title = element_text(face = 2, hjust = 0.5),
panel.border = element_rect(fill = NA, linewidth = 0.2))
GPplot
PLOTS BASED ON MODELS: CALCIFICATION
C_silicate2 <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C") %>%
filter(SGD_number!=0) %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=silicate_umolL,
y=umol.cm2.hr,
color=as.factor(new_colonynumber))) +
geom_point(size=4) +
facet_wrap(~new_colonynumber) +
# scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
theme_classic() +
coord_trans(x ="log") +
geom_smooth(method="lm", formula="y~x + I(x^2)") +
labs(x = expression(Silicate~(mu*mol~L^-1)),
y = expression(Calcification~Rate~(mu*mol~O[2]~cm^2~hr^-1))) +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=20, face="bold"),
axis.title.y=element_text(size=20, face="bold"),
legend.title = element_text(size=18,
face="bold"),
legend.text = element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20,
face="bold"))
C_silicate2 
Patch plots together
###########################
## patch together
###########################
#### remove x-axis labels
p1 <- R_silicate + theme(axis.title.x = element_blank(), axis.text.x = element_blank(),
legend.title = element_blank(), legend.text = element_blank(), legend.position = "none",
axis.title.y=element_text(size=15, face="bold")) +
scale_x_continuous(breaks = c(2, 10, 20))
p2 <- GP_silicate2 + theme(axis.title.x = element_blank(), axis.text.x = element_blank(),
axis.title.y=element_text(size=15, face="bold")) +
scale_x_continuous(breaks = c(2, 10, 20))
p3 <- C_silicate_ECRS_firebrick4 +
theme(legend.title = element_blank(), legend.text = element_blank(), legend.position = "none",
axis.title.x=element_text(size=15, face="bold"),
axis.title.y=element_text(size=15, face="bold")) +
scale_x_continuous(breaks = c(2, 10, 20))
combined_column <- p1/p2/p3 + plot_layout(guides = 'collect')
combined_column +
plot_annotation(tag_levels = 'A') &
theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
plot.tag = element_text(size = 18, face = "bold"))
ggsave(here::here("Outputs", "Results", "Final", "CombinedRates.jpg"),
width=10, height=15)Plots using pred() function
- as of August 2024, still need to figure out how to do properly
- as of May 15th, 2025, ignore because did not use in paper
#####################################
## New attempt: July 29th, 2024
#####################################
### new dfs just for calcification
RespoR_Normalized_Full_withNEC_nutrients_Conly <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr))
### and sort in the right order from lowest silicate to highest
RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted <- RespoR_Normalized_Full_withNEC_nutrients_Conly %>%
arrange(silicate_umolL)
### this is the MIXED MODEL for C
C_bothsites_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL)*site + (1|new_colonynumber),
data=RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted)
###############################################
## the problem is with the predict function
###############################################
### try with Riley's code
new_data <- expand.grid(
silicate_umolL = seq(min(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$silicate_umolL),
max(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$silicate_umolL), length.out = 100),
site = levels(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$site),
new_colonynumber = unique(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$new_colonynumber))
# Predict values using the model
#new_data$predicted <- predict(C_bothsites_silicate_linearmodel, newdata = new_data, re.form = NA)
riley_pred_plot <- new_data %>%
ggplot(aes(x = silicate_umolL,
y = predicted,
color = site)) +
geom_line(size = 1, linetype = "dashed") + # Add dashed trend lines
geom_point(data = RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted,
aes(x = silicate_umolL,
y = umol.cm2.hr,
color = treatment),
size = 3) + # Add original data points
labs(x = "silicate",
y = "rate",
color = "site") +
theme_minimal() +
theme(legend.title = element_text(face = "bold"))
#riley_pred_plot
###############################
## tidyverse code
###############################
# Generate a sequence of silicate_umolL values for prediction
new_data_C_pred <- tibble(
silicate_umolL = seq(min(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$silicate_umolL),
max(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$silicate_umolL), length.out = 100),
site = unique(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$site)[2], # Use the first site as a placeholder
new_colonynumber = unique(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted$new_colonynumber)[5] # Use the first colony as a placeholder
)
# Add predictions to the new data frame
new_data_C_pred <- new_data_C_pred %>%
mutate(predicted = predict(C_bothsites_silicate_linearmodel, newdata = new_data_C_pred, re.form = NA)) # predict without random effects
Pred_plot_C_10 <- RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted %>%
ggplot(aes(x = log(silicate_umolL),
y = umol.cm2.hr,
color = factor(site))) +
geom_point(size = 3) + # Scatter plot of the original data
geom_line(data = new_data_C_pred, aes(y = predicted), color = 'darkgreen', size = 1) + # Regression line
labs(title = "Mixed-Effects Regression Line Plot",
x = "Silicate (umol/L)",
y = "umol/cm^2/hr") + # Labels
theme_minimal()
Pred_plot_C_10
### it's working rn with this code but now I only have one site line ? Comparing Actual vs Predicted Nutrient Values
- these plots are meant to compare the actual env data collected from the dilutions to the predicted values based on calculations from the endmember gw data
ActualvsPredicted_Nutrients_edit_Nitrate <- ActualvsPredicted_Nutrients %>%
pivot_longer(cols = "Actual":"Predicted", names_to = "Value_Type", values_to = "Value") %>%
filter(Parameter=="Nitrate")
ActualvsPredicted_Nutrients_edit_Phosphate <- ActualvsPredicted_Nutrients %>%
pivot_longer(cols = "Actual":"Predicted", names_to = "Value_Type", values_to = "Value") %>%
filter(Parameter=="Phosphate") %>%
filter(Value<1)
#### averaging the values across days per SGD dil to get just one value
ActualvsPredicted_Nutrients_edit_Nitrate_averaged <- ActualvsPredicted_Nutrients_edit_Nitrate %>%
group_by(SGD_number, Site, Parameter, Value_Type) %>%
dplyr::summarize(meanN = mean(Value))
ActualvsPredicted_Nutrients_edit_Phosphate_averaged <- ActualvsPredicted_Nutrients_edit_Phosphate %>%
group_by(SGD_number, Site, Value_Type) %>%
dplyr::summarize(meanN = mean(Value))
#########################
## nitrate
#########################
ActualvsPredicted_Nutrients_Nitrateplot_AVERAGED <- ActualvsPredicted_Nutrients_edit_Nitrate_averaged %>%
ggplot(aes(x=(SGD_number+1),
y=meanN,
color=Value_Type)) +
geom_point(size=4) +
scale_color_manual(values=pnw_palette("Bay", n=2)) +
facet_wrap(~Site) +
# geom_smooth(method="lm", formula="y~log(x)", color="black") +
theme_bw() +
coord_trans(x ="log") +
labs(x = "SGD Percentages",
y = expression(Nitrate+Nitrite~(mu*mol~L^-1)),
# title = "Predicted vs Actual Values by Nitrate",
color = "Nutrient Measurements") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=20),
axis.title.y=element_text(size=20),
legend.title = element_text(size=18),
legend.text = element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
ActualvsPredicted_Nutrients_Nitrateplot_AVERAGED
ggsave(here::here("Outputs", "Results", "Actual_vs_Predicted_Nitrates_byday_AVG.jpg"),
width = 10, height=7)
#########################
## phosphate
#########################
ActualvsPredicted_Nutrients_Phosphateplot_nooutlier_AVERAGED <- ActualvsPredicted_Nutrients_edit_Phosphate_averaged %>%
ggplot(aes(x=(SGD_number+1),
y=meanN,
color=Value_Type)) +
geom_point(size=4) +
scale_color_manual(values=pnw_palette("Bay", n=2)) +
facet_wrap(~Site) +
# geom_smooth(method="lm", formula="y~log(x)", color="black") +
theme_bw() +
coord_trans(x ="log") +
labs(x = "SGD Percentages",
y = expression(Phosphate~(mu*mol~L^-1)),
# title = "Predicted vs Actual Values by Phosphate",
color = "Nutrient Measurements") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=20),
axis.title.y=element_text(size=20),
legend.title = element_text(size=18),
legend.text = element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
ActualvsPredicted_Nutrients_Phosphateplot_nooutlier_AVERAGED
ggsave(here::here("Outputs", "Results", "Actual_vs_Predicted_Phosphate_nooutlier_AVG.jpg"),
width = 10, height=7)
### patch them together
p1.1a <- ActualvsPredicted_Nutrients_Nitrateplot_AVERAGED + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank())
p1.1b <- ActualvsPredicted_Nutrients_Phosphateplot_nooutlier_AVERAGED
actualVpredicted_NutsPlot <- p1.1a/p1.1b +
plot_layout(guides = 'collect') +
plot_annotation(tag_levels = 'A') #+
#theme(axis.title = element_text(face = "bold"),
# plot.tag = element_text(size = 12, face = "bold"))
actualVpredicted_NutsPlot
ggsave(here::here("Outputs", "Results", "Final", "ActualVPredicted_NutsPlot.jpg"),
width = 12, height=8)
#########################################
## look at day by day measurements
## because averages are off for first three dilutions for actual data
#########################################
### NN
NN_byday_C <- ActualvsPredicted_Nutrients_edit_Nitrate %>%
filter(Site=="Cabral") %>%
ggplot(aes(x=(SGD_number+1),
y=Value,
color=Value_Type)) +
geom_point() +
facet_wrap(~Day) +
scale_color_manual(values=pnw_palette("Bay", n=2)) +
coord_trans(x ="log") +
labs(title="NN for Cabral")
NN_byday_C
NN_byday_V <- ActualvsPredicted_Nutrients_edit_Nitrate %>%
filter(Site=="Varari") %>%
ggplot(aes(x=(SGD_number+1),
y=Value,
color=Value_Type)) +
geom_point() +
facet_wrap(~Day) +
scale_color_manual(values=pnw_palette("Bay", n=2)) +
coord_trans(x ="log") +
labs(title="NN for Varari")
NN_byday_V
### Phos
Phos_byday_C <- ActualvsPredicted_Nutrients_edit_Phosphate %>%
filter(Site=="Cabral") %>%
ggplot(aes(x=(SGD_number+1),
y=Value,
color=Value_Type)) +
geom_point() +
facet_wrap(~Day) +
scale_color_manual(values=pnw_palette("Bay", n=2)) +
coord_trans(x ="log") +
labs(title="Phosphate for Cabral")
Phos_byday_C
Phos_byday_V <- ActualvsPredicted_Nutrients_edit_Phosphate %>%
filter(Site=="Varari") %>%
ggplot(aes(x=(SGD_number+1),
y=Value,
color=Value_Type)) +
geom_point() +
facet_wrap(~Day) +
scale_color_manual(values=pnw_palette("Bay", n=2)) +
coord_trans(x ="log") +
labs(title="Phosphate for Varari")
Phos_byday_V
#########################
### Patch
#########################
NN_byday_C + Phos_byday_C + NN_byday_V + Phos_byday_V
ggsave(here::here("Outputs", "Results", "exploring_actualvspred.jpg"),
width = 12, height=7)check GP and C again colored by site and TA by NEC rates
- Results show that TA and NEC relationship is opposite of what we expect … so it is likely that nutrients have a stronger effect on the physiological rates than TA
### GP on x and C on y - colored by site
RespoR_Normalized_Full_withNEC_nutrients2 <- RespoR_Normalized_Full_withNEC_nutrients %>%
select(!c("phosphate_umolL":"TA_initial")) %>%
filter(P_R!="R", P_R!="NP") %>%
pivot_wider(names_from = P_R, values_from = umol.cm2.hr)
GP_bothsites_forNEC <- RespoR_Normalized_Full_withNEC_nutrients2 %>%
filter(sample_ID!="Varari_Col8_Dil6_Light") %>%
# filter(P_R=="C") %>%
filter(!is.infinite(GP)) %>%
ggplot(aes(x=GP,
y=C)) +
geom_point(size=4, aes(color=site)) +
scale_color_manual(values=pnw_palette("Bay", n=2),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~log(x)", color="black") +
geom_hline(yintercept = 0) +
labs(x = "GP",
y = "C",
title = "GP by C both sites") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
GP_bothsites_forNEC
## model
model_GP_C_bothsites <- lmer(GP ~ C*site + (1|new_colonynumber), data = RespoR_Normalized_Full_withNEC_nutrients2 %>%
filter(sample_ID!="Varari_Col8_Dil6_Light") %>%
filter(!is.infinite(GP)))
anova(model_GP_C_bothsites)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## C 0.032975 0.032975 1 112.006 1.9160 0.1691
## site 0.000061 0.000061 1 19.805 0.0035 0.9533
## C:site 0.006054 0.006054 1 112.006 0.3517 0.5543############################
## TA by NEC
###########################
TA_bothsites_forNEC <- RespoR_Normalized_Full_withNEC_nutrients %>%
#filter(sample_ID=="Varari_Col8_Dil6_Light") %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=TA_initial,
y=umol.cm2.hr)) +
geom_point(size=4, aes(color=site)) +
scale_color_manual(values=pnw_palette("Bay", n=2),
guide="none") +
theme_bw() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~log(x)", color="black") +
geom_hline(yintercept = 0) +
labs(x = "TA",
y = "C",
title = "TA by C both sites") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=18),
axis.title.y=element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20))
TA_bothsites_forNEC
For ECRS:
- editing the calcification plots specifically for the ECRS talk
###############################
## Calcification
###############################
C_silicate_ECRS_firebrick4 <- RespoR_Normalized_Full_withNEC_nutrients %>%
mutate(site = dplyr::recode(site,
"Varari" = "Site 2",
"Cabral" = "Site 1")) %>%
filter(P_R=="C") %>%
filter(!is.infinite(umol.cm2.hr)) %>%
ggplot(aes(x=silicate_umolL,
y=umol.cm2.hr)) +
geom_point(size=5, aes(color=site)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4")) +
theme_classic() +
coord_trans(x ="log") +
theme(strip.background = element_rect(fill = "white")) +
geom_hline(yintercept = 0) +
geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R == "C", site == "Varari", !is.infinite(umol.cm2.hr)),
aes(x = silicate_umolL, y = umol.cm2.hr),
method = "lm", formula = y ~ log(x), color = "firebrick4") +
geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R == "C", site == "Cabral", !is.infinite(umol.cm2.hr)),
aes(x = silicate_umolL, y = umol.cm2.hr),
method = "lm", formula = y ~ log(x), color = "darkgoldenrod2", linetype = "dashed") +
# geom_smooth(method="lm", formula = "y~log(x)", color="red3",
# method="lm", formula = "y~log(x)", color="deepskyblue4", type="dashed") +
# geom_smooth(method="lm", formula = "y~log(x)", color="deepskyblue4", type="dashed") +
labs(x = "Logged Silicate (umolL)",
y = "Rate (umol CaCO3 g-1 hr-1)",
title = "Calcification Rates as a Function of Silicate",
color = "Site") +
theme(axis.text.x=element_text(size=18),
axis.text.y=element_text(size=18),
axis.title.x=element_text(size=20, face="bold"),
axis.title.y=element_text(size=20, face="bold"),
legend.title = element_text(size=18,
face="bold"),
legend.text = element_text(size=18),
plot.title=element_text(hjust=0.5, ## centers title
size=20,
face="bold"))
C_silicate_ECRS_firebrick4
ggsave(here::here("Outputs", "Results", "C_silicate_ECRS.jpg"),
width=10, height=7)EXtra:
- not really using and not commented but saving as back up
#RespoR_Normalized_Full_withNEC_nutrients_Cabral_GP <- RespoR_Normalized_Full_withNEC_nutrients %>%
# filter(site=="Cabral") %>%
# filter(P_R=="C") %>%
# filter(!is.infinite(umol.cm2.hr))
#min(RespoR_Normalized_Full_withNEC_nutrients_Cabral_GP$umol.cm2.hr)
#max(RespoR_Normalized_Full_withNEC_nutrients_Cabral_GP$umol.cm2.hr)
#RespoR_Normalized_Full_withNEC_nutrients_Varari_GP <- RespoR_Normalized_Full_withNEC_nutrients %>%
# filter(site=="Varari") %>%
# filter(P_R=="C") %>%
# filter(!is.infinite(umol.cm2.hr))
#min(RespoR_Normalized_Full_withNEC_nutrients_Varari_GP$umol.cm2.hr)
#max(RespoR_Normalized_Full_withNEC_nutrients_Varari_GP$umol.cm2.hr)
### plot
#C_silicate_edit <- RespoR_Normalized_Full_withNEC_nutrients %>%
#filter(P_R=="C") %>%
# filter(!is.infinite(umol.cm2.hr)) %>%
#ggplot(aes(x=silicate_umolL,
# y=umol.cm2.hr)) +
# geom_point(size=4, aes(color=site)) +
# scale_color_manual(values=pnw_palette("Bay", n=2)) +
#theme_classic() +
#coord_trans(x ="log") +
#geom_abline(slope = -0.02219, intercept = 0.15352, color = "black") +
# geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2)
#geom_smooth(data = RespoR_Normalized_Full_withNEC_nutrients %>%
# filter(P_R == "C", site == "Varari", !is.infinite(umol.cm2.hr)),
# aes(x = silicate_umolL, y = umol.cm2.hr),
# method = "lm", formula = y ~ log(x), color = "red3")
# labs(x = "logged Silicate (umolL)",
# y = "Rate (umol CaCO3 g-1 hr-1)",
# title = "Calcification Rates as a Function of Silicate",
# color = "Site") +
# theme(axis.text.x=element_text(size=18),
# axis.text.y=element_text(size=18),
# axis.title.x=element_text(size=20),
# axis.title.y=element_text(size=20),
## legend.title = element_text(size=18),
# legend.text = element_text(size=18),
# plot.title=element_text(hjust=0.5, ## centers title
# size=20))
#C_silicate_edit
##############################################################
##############################################################
###############################
### trying to use marginal and conditional values to plot
###############################
### mixed model
#C_bothsites_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL)*site + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
# filter(P_R=="C") %>%
# filter(!is.infinite(umol.cm2.hr)))
#C_silicate2 <- RespoR_Normalized_Full_withNEC_nutrients %>%
# filter(P_R=="C") %>%
# filter(!is.infinite(umol.cm2.hr)) %>%
# mutate(fit.m = predict(C_bothsites_silicate_linearmodel),
# fit.c = predict(C_bothsites_silicate_linearmodel)) #+
# mutate(resid = resid(C_bothsites_silicate_linearmodel)) ### getting an error here
#C_silicate2.1 <- C_silicate2 %>%
# ggplot(aes(x=silicate_umolL,
# y=umol.cm2.hr)) +
# geom_point(pch = 16, col = "grey") +
#geom_line(aes(y = fit.c), col = 1, size = 2)
##C_silicate2.1
# geom_point(size=4, aes(color=site)) +
# scale_color_manual(values=pnw_palette("Bay", n=2)) +
#theme_bw() +
# coord_trans(x ="log") +
# theme(strip.background = element_rect(fill = "white")) +
# geom_hline(yintercept = 0) +
# geom_line(aes(y = fit.c), col = 1, size = 2) +
# coord_cartesian(ylim = c(-40, 100)) +
#geom_smooth(method="lm", formula = "y~log(x)", color="black") +
# labs(x = "logged Silicate (umolL)",
# y = "Rate (umol O2 cm2 hr-1)",
# title = "Calcification Rates as a Function of Silicate",
# color = "Site") +
# theme(axis.text.x=element_text(size=18),
# axis.text.y=element_text(size=18),
# axis.title.x=element_text(size=18),
# axis.title.y=element_text(size=18),
# plot.title=element_text(hjust=0.5, ## centers title
# size=20))
#C_silicate2.1#
#ggsave(here::here("Outputs", "Results", "C_silicate.jpg"))
###############################
## for Calcification
###############################
#RespoR_Normalized_Full_withNEC_nutrients2<- RespoR_Normalized_Full_withNEC_nutrients %>%
#filter(P_R=="C") %>%
# filter(!is.infinite(umol.cm2.hr)) %>%
# arrange(RespoR_Normalized_Full_withNEC_nutrients2$umol.cm2.hr)
#C_bothsites_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL)*site + (1|new_colonynumber), #data=RespoR_Normalized_Full_withNEC_nutrients2)
### ^^ look for slope and intercept
#summary(C_bothsites_silicate_linearmodel)
### use predict function to get info for confidence intervals - use cbind to create a new df and use that
#predC <- predict(C_bothsites_silicate_linearmodel) # Gets the predicted values from the regression lines in the ANOVA
#graphdata_C <- cbind(RespoR_Normalized_Full_withNEC_nutrients2, predC) # attaches those predictions to the dataset
#graphdata_C <- graphdata_C %>%
#arrange(graphdata_C$silicate_umolL)
# We'll plot the original data as points, and the regression predictions as lines
#C_silicate_pred <- new_data %>%
# ggplot(aes(x=silicate_umolL,
# y=umol.cm2.hr,
# color=site)) +
# geom_point(size=4) +
# geom_line(aes(y=(predC))) +
# geom_line(data = new_data, aes(x = log(silicate_umolL), y = predictions, color = site)) +
# scale_color_manual(values=pnw_palette("Bay", n=2)) +
# theme_classic() +
# coord_trans(x ="log") +
# labs(x = "logged Silicate (umolL)",
# y = "Rate (umol CaCO3 g-1 hr-1)",
# title = "Calcification Rates as a Function of Silicate",
# color = "Site") +
# theme(axis.text.x=element_text(size=18),
# axis.text.y=element_text(size=18),
# axis.title.x=element_text(size=20),
# axis.title.y=element_text(size=20),
# legend.title = element_text(size=18),
# legend.text = element_text(size=18),
# plot.title=element_text(hjust=0.5, ## centers title
# size=20))
#C_silicate_pred
###############################
## for Gross Photosynthesis
###############################
#RespoR_Normalized_Full_withNEC_nutrients3 <- RespoR_Normalized_Full_withNEC_nutrients %>%
# # filter(P_R=="GP") %>%
# filter(!is.infinite(umol.cm2.hr))
#GP_bothsites_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL)*site + (1|new_colonynumber), #data=RespoR_Normalized_Full_withNEC_nutrients3)
### ^^ look for slope and intercept
#summary(GP_bothsites_silicate_linearmodel)
############################
### attemot 1 -- something is off with the pred function, this is for a simpler model, mine is too complex
############################
#Cpred <- predict(C_bothsites_silicate_linearmodel)
#graph_Cpred_data <- cbind(RespoR_Normalized_Full_withNEC_nutrients_Conly_sorted, Cpred)
#Cpred_plot <- ggplot(data=graph_Cpred_data,
# aes(x=log(silicate_umolL),
# y=umol.cm2.hr,
# color=site)) +
# geom_point(size=4) +
# geom_line(aes(y=Cpred)) +
# labs(x="silicate_umolL", y="umol.cm2.hr", fill="site") +
# theme_bw()
#Cpred_plot
#Cpred_plot2 <- ggplot(data=graph_Cpred_data,
# aes(x=log(silicate_umolL),
# y=Cpred,
# color=site)) +
# geom_point(size=4) +
# geom_line(aes(y=Cpred)) +
# labs(x="silicate_umolL", y="Cpred", fill="site") +
# theme_bw()
#Cpred_plot2
###########################
## get basic correlations with p-values
##########################
#res.cor.Varari <- correlate(corr_matrix_Varari, method = "spearman", diagonal = 1)
###########################
## edit correlation matrix and plot
##########################
#res.cor.Varari %>%
#rearrange(method = "MDS", absolute = FALSE) %>%
# shave(upper=TRUE) %>%
# filter(!Salinity = 1) %>%
# rplot() ## make correlologram
###########################
#res.cor.Varari %>%
#filter(!"Salinity" = 1)
# rearrange(method = "MDS", absolute = FALSE) %>%
# shave(upper=TRUE) %>%
##########################Make correlation matrix for nutrients
## data
corr_matrix_envparams <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(P_R=="R") %>%
mutate(Salinity=salinity) %>%
select(!c("salinity", "umol.cm2.hr", "SGD_number", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL")) ## removed ammonia for rn because have less than values and causing NAs
##visualize it
ggcorr(corr_matrix_envparams, method = c("everything", "pearson")) 
###############################
### Dani B's methods #######
cor_mat <- round(cor(corr_matrix_envparams),4)
head(cor_mat)## phosphate_umolL silicate_umolL NN_umolL temp_C new_pH
## phosphate_umolL 1.0000 0.2483 0.2789 0.1626 -0.0260
## silicate_umolL 0.2483 1.0000 0.1296 -0.1146 -0.1927
## NN_umolL 0.2789 0.1296 1.0000 0.1277 -0.0638
## temp_C 0.1626 -0.1146 0.1277 1.0000 0.2110
## new_pH -0.0260 -0.1927 -0.0638 0.2110 1.0000
## TA_initial NA NA NA NA NA
## TA_initial log_silicate log_silicate2 Salinity
## phosphate_umolL NA 0.2378 0.2483 -0.2635
## silicate_umolL NA 0.9590 0.9887 -0.6914
## NN_umolL NA 0.1640 0.1465 0.2955
## temp_C NA -0.1790 -0.1415 -0.1727
## new_pH NA -0.2150 -0.2027 -0.0230
## TA_initial 1 NA NA NA#melt the correlation matrix means it reassembles data frame to be more effective to complete corr matrix
#to long format
melted_cormat <- melt(cor_mat)
head(melted_cormat)## Var1 Var2 value
## 1 phosphate_umolL phosphate_umolL 1.0000
## 2 silicate_umolL phosphate_umolL 0.2483
## 3 NN_umolL phosphate_umolL 0.2789
## 4 temp_C phosphate_umolL 0.1626
## 5 new_pH phosphate_umolL -0.0260
## 6 TA_initial phosphate_umolL NA#visulaize the correlation matrix in general
ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) +
geom_tile()
# Get lower and upper triangle of the correlation matrix
#Note that, a correlation matrix has redundant information. We’ll use the functions below to set half of it to NA
get_lower_tri<-function(cor_mat){
cormat[upper.tri(cor_mat)] <- NA
return(cor_mat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cor_mat){
cor_mat[lower.tri(cor_mat)]<- NA
return(cor_mat)
}
#apply upper tri calculation to graphc
upper_tri <- get_upper_tri(cor_mat)
upper_tri## phosphate_umolL silicate_umolL NN_umolL temp_C new_pH
## phosphate_umolL 1 0.2483 0.2789 0.1626 -0.0260
## silicate_umolL NA 1.0000 0.1296 -0.1146 -0.1927
## NN_umolL NA NA 1.0000 0.1277 -0.0638
## temp_C NA NA NA 1.0000 0.2110
## new_pH NA NA NA NA 1.0000
## TA_initial NA NA NA NA NA
## log_silicate NA NA NA NA NA
## log_silicate2 NA NA NA NA NA
## Salinity NA NA NA NA NA
## TA_initial log_silicate log_silicate2 Salinity
## phosphate_umolL NA 0.2378 0.2483 -0.2635
## silicate_umolL NA 0.9590 0.9887 -0.6914
## NN_umolL NA 0.1640 0.1465 0.2955
## temp_C NA -0.1790 -0.1415 -0.1727
## new_pH NA -0.2150 -0.2027 -0.0230
## TA_initial 1 NA NA NA
## log_silicate NA 1.0000 0.9899 -0.5874
## log_silicate2 NA NA 1.0000 -0.6489
## Salinity NA NA NA 1.0000#melt the correlation matrix
#melt the correlation data and drop the rows with NA values
melted_cormat <- melt(upper_tri, na.rm = TRUE)
#heatmap of correlation matrix
#negative correlations are in purple color and positive correlations in red
#scale_fill_gradient2 is used with the argument limit = c(-1,1) as correlation coefficients range from -1 to 1
#coord_fixed() : this function ensures that one unit on the x-axis is the same length as one unit on the y-axis
ggplot(data = melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "midnightblue", high = "firebrick4", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
melted_cormat <- melted_cormat %>% mutate_at(vars(starts_with("value")), funs(round(., 2)))
# Create a ggheatmap with basic characteristics, etc.
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "firebrick3", mid = "white", high = "dodgerblue3",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
# Print the heatmap
print(ggheatmap)
#add correlation coefficients to the heatmap
#geom_text() to add the correlation coefficients on the graph
#guides() to change the position of the legend title
#if else statement in melted data frame to quotes of black and white to adjust text color
corr_matrix_SGD_params <- ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 6) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(size = 18, face="bold", color="black"),
axis.text.y = element_text(size = 18, face="bold", color="black"),
panel.grid.major = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.ticks = element_blank(),
legend.justification = c(0.8, 0),
legend.title = element_text(size = 18, face="bold", color="black"),
legend.text = element_text(size = 20, face="bold", color="black"),
legend.position = c(0.48, 0.75),
legend.direction = "horizontal") +
guides(fill = guide_colorbar(barwidth = 12, barheight = 2,
title.position = "top", title.hjust = 0.5, title.vjust = 1.0))
corr_matrix_SGD_params
Figure 2:
############################################
### correlation matrix for Varari
############################################
## data
corr_matrix_envparams_Varari <- RespoR_Normalized_Full_withNEC_nutrients %>%
mutate(N_P = NN_umolL/phosphate_umolL) %>%
filter(P_R=="R") %>%
mutate(Salinity=salinity,
TA = TA_initial,
pH = new_pH,
Nitrate_Nitrite_umolL = NN_umolL,
Phosphate_umolL = phosphate_umolL,
Silicate_umolL = silicate_umolL,
Ammonia_umolL = ammonia_umolL) %>%
filter(site=="Varari") %>%
select(!c("salinity", "silicate_umolL", "TA_initial", "Ammonia_umolL", "new_pH", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL"))
##visualize it
ggcorr(corr_matrix_envparams_Varari, method = c("everything", "pearson")) 
###############################
cor_mat <- round(cor(corr_matrix_envparams_Varari),4)
#head(cor_mat)
###############################
# Set the diagonal elements to NA
#diag(cor_mat) <- NA
###############################
###############################
# Melt the correlation matrix
melted_cormat <- melt(cor_mat, na.rm = TRUE)
###############################
#melt the correlation matrix means it reassembles data frame to be more effective to complete corr matrix
#to long format
#melted_cormat <- melt(cor_mat)
#head(melted_cormat)
#visualize the correlation matrix in general
#ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) +
# geom_tile()
# Get lower and upper triangle of the correlation matrix
#Note that, a correlation matrix has redundant information. We’ll use the functions below to set half of it to NA
get_lower_tri<-function(cor_mat){
cormat[upper.tri(cor_mat)] <- NA
return(cor_mat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cor_mat){
cor_mat[lower.tri(cor_mat)]<- NA
return(cor_mat)
}
#apply upper tri calculation to graphc
upper_tri <- get_upper_tri(cor_mat)
#upper_tri
#melt the correlation matrix
#melt the correlation data and drop the rows with NA values
melted_cormat <- melt(upper_tri, na.rm = TRUE)
#heatmap of correlation matrix
#negative correlations are in purple color and positive correlations in red
#scale_fill_gradient2 is used with the argument limit = c(-1,1) as correlation coefficients range from -1 to 1
#coord_fixed() : this function ensures that one unit on the x-axis is the same length as one unit on the y-axis
ggplot(data = melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "midnightblue", high = "firebrick4", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson Correlation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
melted_cormat <- melted_cormat %>% mutate_at(vars(starts_with("value")), funs(round(., 2)))
# Create a ggheatmap with basic characteristics, etc.
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "firebrick3", mid = "white", high = "dodgerblue3",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson Correlation") +
theme_minimal()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 10, hjust = 1))+
coord_fixed()
diag(melted_cormat) <- NA
# Print the heatmap
#print(ggheatmap)
#add correlation coefficients to the heatmap
#geom_text() to add the correlation coefficients on the graph
#guides() to change the position of the legend title
#if else statement in melted data frame to quotes of black and white to adjust text color
corr_matrix_SGD_params_Varari <- ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(size = 12, face="bold", color="black"),
axis.text.y = element_text(size = 12, face="bold", color="black"),
panel.grid.major = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.ticks = element_blank(),
legend.justification = c(0.8, 0),
legend.title = element_text(size = 10, color="black"),
legend.text = element_text(size = 10, color="black"),
legend.position = c(0.48, 0.75),
legend.direction = "horizontal",
title= element_text(size = 12, hjust=0.5, face="bold", color="black")) +
labs(title="Pearson Correlation Matrix for Environmental Parameters at Varari") +
guides(fill = guide_colorbar(barwidth = 8, barheight = 1,
title.position = "top", title.hjust = 0.5, title.vjust = 1.0))
corr_matrix_SGD_params_Varari
ggsave(here::here("Outputs", "Results", "Varari_correlation_matrix.jpg"),
width=10, height=7)Figure 2a:
############################################
### correlation matrix for Cabral
############################################
## data
corr_matrix_envparams_Cabral <- RespoR_Normalized_Full_withNEC_nutrients %>%
mutate(N_P = NN_umolL/phosphate_umolL) %>%
filter(P_R=="R") %>%
mutate(Salinity=salinity,
TA = TA_initial,
pH = new_pH,
Nitrate_Nitrite_umolL = NN_umolL,
Phosphate_umolL = phosphate_umolL,
Silicate_umolL = silicate_umolL,
Ammonia_umolL = ammonia_umolL) %>%
filter(site=="Cabral") %>%
select(!c("salinity", "silicate_umolL", "TA_initial", "new_pH", "phosphate_umolL", "NN_umolL", "umol.cm2.hr", "SGD_number", "temp_C", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL", "Ammonia_umolL"))
##visualize it
ggcorr(corr_matrix_envparams_Cabral, method = c("everything", "pearson")) 
###############################
cor_mat <- round(cor(corr_matrix_envparams_Cabral),4)
#head(cor_mat)
#melt the correlation matrix means it reassembles data frame to be more effective to complete corr matrix
#to long format
melted_cormat <- melt(cor_mat)
#head(melted_cormat)
#visualize the correlation matrix in general
ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) +
geom_tile()
# Get lower and upper triangle of the correlation matrix
#Note that, a correlation matrix has redundant information. We’ll use the functions below to set half of it to NA
get_lower_tri<-function(cor_mat){
cormat[upper.tri(cor_mat)] <- NA
return(cor_mat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cor_mat){
cor_mat[lower.tri(cor_mat)]<- NA
return(cor_mat)
}
#apply upper tri calculation to graphc
upper_tri <- get_upper_tri(cor_mat)
#upper_tri
#melt the correlation matrix
#melt the correlation data and drop the rows with NA values
melted_cormat <- melt(upper_tri, na.rm = TRUE)
#heatmap of correlation matrix
#negative correlations are in purple color and positive correlations in red
#scale_fill_gradient2 is used with the argument limit = c(-1,1) as correlation coefficients range from -1 to 1
#coord_fixed() : this function ensures that one unit on the x-axis is the same length as one unit on the y-axis
ggplot(data = melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "midnightblue", high = "firebrick4", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson Correlation") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
melted_cormat <- melted_cormat %>% mutate_at(vars(starts_with("value")), funs(round(., 2)))
# Create a ggheatmap with basic characteristics, etc.
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "firebrick3", mid = "white", high = "dodgerblue3",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson Correlation") +
theme_minimal()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 10, hjust = 1))+
coord_fixed()
# Print the heatmap
#print(ggheatmap)
#add correlation coefficients to the heatmap
#geom_text() to add the correlation coefficients on the graph
#guides() to change the position of the legend title
#if else statement in melted data frame to quotes of black and white to adjust text color
corr_matrix_SGD_params_Cabral <- ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(size = 12, face="bold", color="black"),
axis.text.y = element_text(size = 12, face="bold", color="black"),
panel.grid.major = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.ticks = element_blank(),
legend.justification = c(0.8, 0),
legend.title = element_text(size = 10, color="black"),
legend.text = element_text(size = 10, color="black"),
legend.position = c(0.48, 0.75),
legend.direction = "horizontal",
title= element_text(size = 12, hjust=0.5, face="bold", color="black")) +
labs(title="Pearson Correlation Matrix for Environmental Parameters at Cabral") +
guides(fill = guide_colorbar(barwidth = 8, barheight = 1,
title.position = "top", title.hjust = 0.5, title.vjust = 1.0))
corr_matrix_SGD_params_Cabral
ggsave(here::here("Outputs", "Results", "Cabral_correlation_matrix.jpg"),
width=10, height=7)Results of back calculating
Essentially found that the predicted values were too low to be useful for
take average of Varari Days 2-4 to fill in for ambient day 1
- Did not end up doing, May 15th, 2025
C_only <- RespoR_Normalized_Full_withNEC_nutrients %>%
filter(!is.infinite(umol.cm2.hr)) %>%
filter(P_R=="C") %>%
mutate(new_colonynumber=as.numeric(new_colonynumber))
-18.427160## [1] -18.42716-57.911680## [1] -57.91168-63.030418## [1] -63.03042(-18.427160+-57.911680+-63.030418)/3## [1] -46.45642## REPLACE 146 FROM AMBIETN DAY 1 WITH -46.45642 WHICH IS AVERAGE OF ALL 3 OTHER DAYS FOR DELTA TA AT SGD NUMBER 1, AMB SW